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1 .  INTRODUCTION 


1.1  Identification  of  the  Problem 

The  U.S.  Army  currently  has  the  need  for  the  capability  of  rapidly  modeling  and 
evaluating  the  multispectral  obscuration  capacity  of  a  smoke  cloud  in  space  and 
time.  Such  a  multispectral  smoke  obscuration  model  requires  that  the  fate  of  a 
number  of  particle  size  regimes  be  accounted  for  explicitly,'  as  the  scattering 
cross  section  of  a  smoke  particle  at  a  specific  incident  wavelength  depends 
strongly  on  the  particle's  diameter.  The  further  constraints  that  these 
simulations  include  realistic  three-dimensional  mesoscale  structures  of  the  flow 
and  turbulence  fields,  as  well  as  resolve  the  time  dependence  of  the  obscuring 
cloud's  behavior  on  "relevant"  time  scales,  suggest  a  model  of  substantial 
sophistication  and  computational  intensity.  Nevertheless,  a  further  objective 
is  to  perform  such  modeling  using  the  rather  modest  resources  of  a  modern 
personal  computer  (PC) . 

1.2  Development  of  an  Approach 

A  typical  mesoscale  wind  field  model  may  cover  a  lateral  domain  of  20  to  200  km 
with  a  typical  resolution  of  from  500  m  to  5  km  and  a  vertical  domain  of  up  to 
5  km  with  resolutions  ranging  from  several  meters  near  the  surface  to  a  kilometer 
or  more  near  the  top  of  the  domain.  Such  spatial  resolution  is  often  far  greater 
than  the  dimensions  of  the  smoke  cloud  being  modeled  and,  of  course,  only 
reflects  those  spatial  structures  in  the  flow  having  wavelengths  greater  than  the 
resolution  of  the  mesh.  In  addition,  the  temporal  resolution  of  these  wind 
fields  is  typically  in  the  range  of  fractions  of  an  hour  (for  example,  10 
minutes)  to  one  to  several  hours.  Such  time  resolutions  may  be  adequate  to 
resolve  the  evolution  of  the  mean  flow  but  are  clearly  too  coarse  to  resolve  the 
influence  of  intermittent  moderate  size  (for  example,  of  the  size  order  of  the 
smoke  cloud)  sub -grid- scale  eddies  on  the  smoke  cloud. 

Two  possible  modeling  approaches  immediately  come  to  mind- -large  eddy  simulation 
(LES)  and  Monte  Carlo  Lagrangian  particle  dispersion.  The  high  space- time 
resolution  of  an  LES  model  could  be  used  to  drive  a  particle  or  Eulerian  grid 
model;  however,  the  required  computational  intensity  of  LES  models  calls  for  a 
supercomputer  rather  than  a  PC.  Alternatively,  one  could  use  the  wind 
field/boundary  layer  model's  specification  of  the  local  turbulence  fields  to 
drive  a  Monte  Carlo  Lagrangian  particle  model.  Such  an  approach  generally 
considers  each  particle  (or  particle  pair)  to  move  independently  and,  therefore, 
leads  to  ensemble  average  concentration  statistics,  which  for  stationary 
turbulence  can  be  considered  equivalent  to  time -averaged  concentrations.  This 
relatively  computationally  inexpensive  modeling  approach  is  excellent  for  many 
applications  but  could  yield  misleading  results  in  an  environment  where  the  issue 
of  seeing  and  being  seen  may  involve  brief  periods  of  time  rather  than  tirae- 
averaged  conditions. 

A  compromise  approach  might  involve  the  use  of  "canned"  sub-grid  scale  flow 
fields  developed  in  advance  with  the  aid  of  LES  model  simulations,  combined  with 
a  Lagrangian  particle  model  designed  to  let  particles  follow  these  large,  yet 
sub-grid,  scale  motions  augmented  by  a  more  traditional  Monte  Carlo  trajectory 
component  to  account  for  smaller,  unresolved  scale  eddies.  Such  an  approach  has 
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also  been  suggested  by  Sakai  et  al.  (1991),  and  makes  use  of  the  rapidly  evolving 
discipline  of  kinematic  simulation.  This  methodology  will  serve  as  the  basis  for 
moving  particles  within  the  proposed  model. 

Of  course,  it  is  not  enough  to  simply  move  the  particles  correctly;  we  must  also 
account  for  their  size-dependent  removal  rates  and  determine  cost  effective  ways 
to  assess  local  particle  concentrations  and  path-averaged  measures.  Both  of 
these  issues  can  have  very  large  computing  cost  implications  for  a  model.  For 
example,  particles  which  are  effective  scatterers  in  the  far  infrared  (IR)  (for 
example,  particle  diameter  d  »  wavelength  1  «  100  /i)  can  settle  out  with 
deposition  velocities  of  order  10  cm/s  whereas  optimal  ultraviolet  (UV) 
scatterers  can  have  deposition  velocities  ranging  from  10  ^  to  10"*^  cm/s  but  all 
negligible  relative  to  turbulent  velocities  of  order  5-50  cm/s.  Thus,  it  may  be 
possible  to  simulate  all  particle  sizes  relevant  to  the  UV  transmission  problem 
with  a  single  mathematical  particle  whereas  the  IR  transmission  problem  may 
require  separate  computation  of  trajectories  of  particles  in  the  various  size 
ranges.  A  truly  comprehensive  multispectral  model  might  therefore  require  a 
choice  (that  is,  user  or  model  selected)  of  wavelength  operational  modes  to  avoid 
excessive  calculation  (for  example,  one  mathematical  particle  for  each  physical 
particle  size  range)  when  a  simpler  calculation  is  adequate. 

Choice  of  an  appropriate  particle  counting  methodology  for  concentration 
evaluation  is  also  essential  in  balancing  the  cost  of  increased  number  of 
particles  with  the  statistical  uncertainties  associated  with  too  few  particles. 
The  proposed  model  will  incorporate  a  modified  "kernel  estimator"  that  diffuses 
either  the  receptor  location  over  a  time- independent  length  scale  or  the  point 
particles  over  a  time -dependent  scale  corresponding  to  the  root  mean  square  (rms) 
separation  generated  by  the  unresolved  fine  scale  eddies. 

1.3  Overall  Model  Design 

Figure  1  displays  the  basic  modules  making  up  the  envisioned  multispectral  smoke 
dispersion  model  and  the  relationship  of  this  model  to  external  data  files  and 
the  typical  user.  Perhaps  the  most  important  input  data  involves  the  three- 
dimensional  gridded  mesoscale  wind  fields  and  accompanying  two-dimensional  fields 
specifying  boundary  layer  quantities  of  interest  (for  example,  mixed  layer 
depths,  friction  velocities ,  heat  fluxes,  convective  velocity  scale).  These  data 
are  provided  by  a  wind  field/PBL  model  external  to  and  separate  from  the 
dispersion  model.  While  the  current  prototype  accesses  data  produced  by  the 
CALMET  meteorological  model  (Scire  et  al.,  1990),  there  are  few  obstacles  in 
substituting  this  model  with  one  potentially  more  suitable  to  the  Army's  needs. 

The  principal  user  interface  via  a  user  control  file  will  ensure  relatively 
uncomplicated,  straightforward  use  of  the  model.  The  user  will  need  to  specify 
the  parameters  associated  with  the  proposed  smoke  release  (for  example,  release 
locations,  release  quantities,  and  proposed  sampling  locations  or  transmission 
paths  of  interest).  Supplementary  detailed  information,  such  as  specification 
of  the  initial  spectral  profile  of  the  smoke  release,  will  be  included  in 
additional  data  files  prepared  well  in  advance  of  actual  applications. 
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Four-Dlaenslonal  Mesoscale 
Multi-Spectral  Smoke 
Dispersion  Model 


Figure  1.  Overall  model  design. 

1.4  Overview  of  the  Technical  Report 

Section  2  of  this  report  reviews  the  current  status  of  mesoscale  models  of  the 
flow  and  turbulence  fields  and  describes  how  this  information  is  incorporated 
into  current  formulations  of  Lagrangian  particle  models  or  the  hybrid  kinematic 
simulation/Langevin  equation  model  envisioned  here.  Section  3  then  deals  with 
those  aspects  of  the  problem  which  involve  the  multispectral  smoke  aspects  of  the 
problem  including  the  dry  deposition  module,  adaptation  of  the  dry  deposition 
velocities  to  evaluate  Lagrangian  particle  capture  probabilities,  and  computation 
of  the  point  and  path  averaged  quantities  of  interest. 

A  number  of  Che  basic  components  needed  for  such  a  dispersion  model  have  been 
coded  and  tested  during  the  course  of  this  project.  These  tests  and  preliminary 
evaluations  are  described  in  section  4. 

Finally,  the  conclusions  and  recommendations  resulting  from  this  Phase  I  research 
and  development  effort  are  presented  in  section  5.  A  complete  list  of  references 
is  also  provided. 
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2.  MESOSCALE  FLOW,  TURBULENCE,  AND  LAGRANGIAN  MODELING  OF  PARTICLE  MOVEMENT 


2.1  Overview  of  Mean  Flow  Models  Incorporating  Complex  Terrain 

A  variety  of  wind  field  generators  is  currently  available,  ranging  from  analytic 
trajectory  generators  for  specific  terrain  features,  to  objective  analysis 
models,  to  diagnostic  generators,  and  ultimately,  to  prognostic  wind  field 
generators  of  varying  complexity.  In  the  following  subsections,  we  describe 
several  of  the  more  readily- available  modeling  approaches  for  each  type  of  wind 
field  generator. 

2.1.1  Objective  analysis  models 

Wind  field  models  in  this  group  use  objective  methods  to  create  an  entire  grid 
of  three-dimensional  wind  vectors  from  an  irregular  spatial  distribution  of 
measured  two-dimensional  vectors.  This  is  typically  done  in  two  steps. 
Initially,  wind  vectors  in  the  horizonta’  are  specified  at  all  points  in  the  grid 
by  interpolating  among  the  measured  wind  vectors.  Once  this  is  completed,  the 
second  step  makes  use  of  procedures  for  modifying  the  vectors  that  were 
interpolated,  with  the  aim  of  minimizing  the  divergence  in  the  flow  field.  The 
influence  of  terrain  is  represented  in  boundary  conditions  that  do  not  .allow  the 
flow  to  penetrate  terrain  features  large  enough  to  be  resolved  by  the  spacing  of 
the  grid  in  the  horizontal.  As  such,  these  models  create  mass -consistent  wind 
fields  representative  of  the  measured  data,  given  the  constraints  imposed  by  the 
topography . 

In  a  recent  review  article,  Hunt  et  al.  (1990),  report  that  this  modeling 
approach  produces  satisfactory  results  for  scales  of  5  km  to  100  km,  especially 
when  wind  data  are  obtained  over  length  scales  consistent  with  changes  in  surface 
topography.  In  fact,  if  the  observed  wind  data  sufficiently  resolve  important 
features  of  the  flow,  then  the  objective  analysis  method  is  likely  to  provide  the 
most  realistic,  and  the  most  cost  effective  means  of  estimating  the  wind  field 
throughout  the  modeling  region. 

Two  well-known  objective  analysis  procedures  are  described  by  Sherman  (1978)  and 
Goodin  et  al.  (1980).  The  first  procedure,  known  as  MATHEW  (Mass-Adjusted 
Three-Dimensional  Wind  Field),  develops  a  three-dimensional  mass -consistent  wind 
field  by  means  of  a  variational  calculus  approach.  The  coordinate  system  is  a 
fixed  grid  which  does  not  follow  terrain.  Terrain  is  represented  in  the  form  of 
obstacle  cells  which  impose  vertical  boundaries  through  which  no  flow  is  allowed. 
In  step  1  of  the  procedure,  interpolation  of  winds  measured  near  the  surface  use 
1/r^  weighting,  while  a  synoptic  analysis  must  be  used  to  provide  the  wind  field 
aloft.  Goodin  et  al.  (1980)  point  out  chat  their  method  incorporates  a 
terrain-following  coordinate  system  with  variable  spacing  in  the  vertical. 
Interpolation  at  the  surface  makes  use  of  1/r^  weighting,  whereas  winds, 
temperatures,  and  mixing  heights  in  the  upper  layers  make  use  of  1/r  weighting. 
Large-scale  terrain  imposes  barriers  to  the  interpolation.  The  influence  of 
sraaller-scale  terrain  (less  than  one  grid  cell  length)  is  included  by  solving 
Poisson's  equation  for  a  forcing  function  that  is  based  on  the  thickness  of  the 
grid-layer,  and  terrain  gradients.  In  the  example  application,  Goodin  et  al . 
(1980)  indicate  that  their  method  is  faster  than  MATHEW,  but  results  in  larger 
horizontal  divergence,  and  therefore  larger  estimates  of  vertical  velocity. 


2.1.2  Diagnostic  models 


Diagnostic  models  are  similar  to  the  objective  analysis  models  in  that  they  rely 
heavily  on  observed  winds  in  order  to  obtain  satisfactory  estimates  of  the  wind 
field  throughout  the  modeling  domain.  However,  in  addition  to  the  requirement 
that  the  modeled  wind  field  be  mass-consistent,  they  typically  introduce  other 
information  about  the  likely  structure  of  the  wind  field  through  use  of 
diagnostic  equations. 

Endlich  et  al.  (1982)  describe  a  diagnostic  model  that  makes  use  of  the  MATHEW 
method  for  reducing  divergence.  They  introduce  geostrophic  winds  above  the 
boundai^  layer,  and  employ  a  logarithmic  vertical  profile  to  interpolate  winds 
between  the  geostrophic  wind  at  the  top  of  the  boundary  layer  and  the  surface 
winds.  They  cast  the  continuity  equation  in  a  terrain- following  coordinate 
system,  and  formulate  the  Lagrange  multiplier  relations.  By  neglecting  some 
terms  that  may  be  significant  in  areas  of  steep  terrain,  the  system  becomes 
analogous  to  Sherman's  (1‘  ^8),  so  the  MATHEW  method  may  be  used.  This  approach 
essentially  incorporates  the  influence  of  terrain  on  the  flow  by  imposing  the 
depth  of  the  boundary  layer  over  the  terrain,  and  forcing  the  flow  within  the 
boundary  layer  to  satisfy  the  continuity  equation. 

A  second  type  of  diagnostic  wind  field  model  is  described  by  Douglas  and  Kessler 
(1988).  Although  it  uses  the  Goodin  et  al.  (1980)  method  for  minimizing 
divergence,  it  first  includes  diagnostic  equations  for  estimating  kinematic 
effects  of  the  terrain,  slope  flows,  and  blocking  effects  of  terrain  in  the 
presence  of  a  stable  density  stratification.  A  mean  wind  for  the  region, 
obtained  from  upper- air  information,  is  used  to  compute  the  vertical  velocity 
associated  with  terrain  features;  in  a  terrain- following  coordinate  system  (Liu 
and  Yocke ,  1980).  The  divergence -minimization  procedure  is  then  applied  to 
modify  the  initial  mean  wind  field.  Slope  flows  estimated  by  procedures 
suggested  by  Allwine  and  Whiteman  (1988)  are  then  added  to  the  field.  Finally, 
if  the  local  Froude  nvunber  at  a  grid  point  is  less  than  a  critical  value,  the 
flow  at  the  grid  point  must  not  have  an  uphill  component.  If  it  shou?d  have  an 
uphill  value,  the  vector  is  turned  so  that  it  is  tangential  to  the  height- contour 
of  the  terrain.  This  stipulation  simulates  the  largely  horizontal  flow  observed 
below  a  "critical  dividing- streamline . " 

Once  this  flow- field  is  complete,  it  embodies  several  terrain- induced  features 
that  may  or  may  not  be  represented  in  the  observed  wind  field.  The  second  step 
in  this  procedure  introduces  the  observed  wind  data.  Values  of  the  windspeed  at 
grid  points  near  monitoring  stations  are  highly  influenced  by  the  obrerved  winds. 
However,  grid  points  that  are  relatively  far  from  monitoring  stations  retain  the 
characteristics  imposed  by  the  diagnostic  relations.  After  smoothing  to  reduce 
discontinuities,  and  adjustments  to  limit  the  vertical  velocity  at  the  top  of  the 
domain,  the  divergence  minimization  procedure  of  Goodin  et  al .  (1980)  is  applied 
to  obtain  the  mass-consistent  wind  field.  The  recently  developed  CALMET  model 
(Scire  et  al.  ,  1990)  incorporates  the  Douglas  and  Kessler  (1988)  diagnostic  model 
together  with  a  complete  model  for  the  PBL  and  micrometeorological  quantities. 

Several  recently  developed  diagnostic  wind  field  models  also  incorporate  the 
concept  of  a  dividing-streamline  or  Froude  number  in  their  formulation  in  order 
to  better  represent  blocking  or  channeling  effects  of  terrain  during  periods  with 
stably-stratified  density  profiles  in  the  vertical.  MUATMOS  (Rosa  et  al .  ,  1988) 
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is  an  adapcacion  of  ATMOSl  (Davis  et  al.,  1984)  which  is  based  on  a  variational 
calculus  method  similar  in  concept  to  MATHEW.  However,  NUATMOS  makes  use  of 
terrain- following  coordinates,  and  includes  a  method  for  prescribing  the 
parameter  a,  which  controls  the  adjustments  made  to  the  vertical  winds  relative 
to  the  horizontal  winds,  by  means  of  the  local  Froude  number.  Ross  et  al.  (1988) 
demonstrate  the  ability  of  NUATMOS  to  characterize  potential  flow  solutions  for 
flow  over  simple  isolated  features  where  a  -  1,  and  discuss  the  success  of 
several  approaches  to  formulating  the  a  parameter  in  terms  of  the  Froude  number 
by  comparing  results  with  tow-tank  results  obtained  by  Hunt  and  Snyder  (1980). 
When  a  -  0,  the  divergence  in  the  flow-field  is  minimized  by  altering  only  the 
horizontal  winds.  By  equating  this  limit  with  the  low  Froude  number  regime,  a 
largely  two-dimensional  flow  is  obtained.  Because  a  local  Froude  number  is  used 
for  each  layer  in  the  model,  the  method  can  result  in  channeled,  two-dimensional 
flow  in  the  lower  layers,  and  nearly  potential  flow  in  the  upper  layers. 

The  WOCSS  (Winds  on  Critical  Streamline  Surfaces)  method  developed  by  Ludwig  and 
Endlich  (1988)  takes  a  substantially  different  approach  to  incorporating 
channeling  effects  that  result  from  stratification.  This  approach  constrains  the 
analysis  by  defining  coordinate  surfaces  that  represent  surfaces  through  which 
no  flow  is  allowed.  These  coordinate  surfaces  may  intersect  the  terrain,  and  are 
determined  objectively  by  essentially  allowing  a  vertical  deflection  caused  by 
terrain  to  be  no  greater  than  the  locar  length  scale,  U/N.  This  means  that 
coordinate  surfaces  will  intersect  terrain  whenever  the  flow  cannot  pass  over  the 
terrain.  Winds  at  grid  points  that  end  up  "inside"  the  terrain  are  set  to  zero, 
thereby  forcing  the  flow  outside  of  the  terrain  to  travel  around  the  terrain 
because  within  each  layer  no  divergenc.s  is  allowed.  This  method  essentially 
turns  the  three-dimensional  flow  problem  into  a  number  of  two-dimensional  flow 
problems.  Mass  fluxes  are  adjusted  to  satisfy  the  continuity  equation  by  means 
of  the  procedure  described  by  Endlich  (1967),  except  that  vorticity  constraints 
are  not  imposed. 

Ludwig  et  al .  (1991)  point  out  that  this  WOCSS  method  does  not  provide  sufficient 
guidance  for  defining  coordinate  surfaces  in  the  absence  of  stable  stratifica¬ 
tion.  Unlike  the  methods  based  on  the  MATHEW  approach,  it  does  not  provide  flow 
fields  similar  to  potential  flow  in  the  absence  of  stratification.  They  suggest 
that  further  development  of  the  method  might  include  combining  it  with  a  two- 
dimensional  prognostic  model  to  produce  realistic  vertical  profiles  of  winds  and 
temperatures  which  can  be  treated  as  "obseirvations . "  Such  an  approach  would 
introduce  diabatic  effects  such  as  drainage  flows  and  coastal  circulations. 

2.1.3  Analytic  flow  models 

When  complex  flow  fields  must  be  resolved  on  scales  of  order  1  km  or  less,  with 
few  wind  measurements  available.  Hunt  et  al.  (1990)  suggest  that  interpolative 
approaches  like  objective  or  diagnostic  procedures  cannot  be  relied  upon.  In 
their  place,  they  advocate  the  use  of  what  we  will  call  analytic  flow  modelo. 
These  approaches  typically  use  scale  analyses  to  simplify  the  equations  of  motion 
within  regions  that  "re  important  to  the  dynamics  governing  the  flow.  These 
regions  are  then  combined  using  the  methods  of  matched  as3nnptotics ;  that  is, 
solutions  for  different  regions  must  match  at  the  interface  and  .lave  similar 
characteristics  at  infinity. 
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In  a  recent  review  paper,  Carruchers  and  Hunt  (1990)  describe  various  flow 
regimes  associated  with  terrain,  identifying  major  dimensionless  parameters  chat 
are  able  to  classify,  or  order,  these  regimes  (after  Hunt  and  Richards,  1984; 
Snyder,  1985).  They  consider  stratification  effects  (neutral  -  to  -  stable) , 
rotation  effects,  and  changes  in  roughness.  However,  they  do  not  consider 
katabatic  effects  such  as  drainage  flows  at  night  or  flows  resulting  from 
differential  heating  during  the  day. 

The  modeling  approach  to  these  flows  is  based  upon  linear  analysis.  That  is,  the 
presence  of  terrain  merely  perturbs  the  flow,  so  that  perturbation  velocities  are 
small  compared  to  the  mean  flow.  Upon  making  this  assumption,  linearized 
versions  of  the  governing  equations  can  be  solved  analytically.  The  most 
complete  approach  of  this  type  is  the  linear  three-dimensional  theory  of  Hunt  et 
al.  (1988).  Carruthers  and  Hunt  (1990)  point  out  chat  this  is  based  on  earlier 
works  of  Jackson  and  Hunt  (1975)  and  Townsend  (1966).  Two  regions  are  defined- - 
an  inner  region  in  which  changes  in  shear  stress  are  important,  and  an  outer 
region  in  which  the  flow  can  be  treated  as  inviscid.  The  depth  of  the  inner 
region  (i)  is  given  by  the  implicit  relation 

{  In  ((/zj  =  2  k^L  a) 

where  is  the  surface  roughness  length,  k  is  von  Karman's  constant,  and  L  is 
the  length  scale  of  the  hill  (see  figure  2).  Throughout  the  inner  region,  the 
turbulence  is  assumed  to  be  in  local  equilibrium.  A  thin  surface  layer  allows 
the  solution  to  be  matched  to  the  surface  boundary  conditions.  The  outer  region 
is  composed  of  two  layers;  a  middle  layer  and  an  outer  layer.  Shear  in  the 
approach  flow  dominates  the  solution  for  the  middle  layer.  The  solution  in 
the  outer  layer  characterizes  the  response  of  the  fluid  to  the  underlying 
perturbations.  The  perturbations  decay  with  height  if  the  stratification  is  very 
weak  (potential  flow  solutions  are  obtained),  and  the  perturbations  produce  waves 
for  moderate  to  strong  stratification. 

The  solutions  for  the  perturbation  velocities  in  the  two  regions  are  obtained  for 
arbitrary  terrain  by  representing  the  terrain  by  its  Fourier  transform.  However, 
the  basic  assumption  that  the  perturbations  are  small  compared  to  the  incident 
flow  requires  that  the  terrain  be  characterized  by  small  slopes,  so  that  the  use 
of  "arbitrary"  terrain  must  be  qualified  accordingly.  The  most  general  numerical 
model  that  is  based  on  these  methods  is  FLOWSTAR  (Carruthers  et  al .  ,  1988). 
Carruthers  and  Hunt  (1990)  point  out  that  the  use  of  Fast  Fourier  Transforms 
(FFT's)  allows  rapid  computations  to  be  made,  and  since  the  method  is  analytic, 
the  flow  field  may  be  found  at  precise  locations,  rather  than  within  a  grid. 
They  also  point  out  that  perturbations  due  to  changes  in  roughness  and  terrain 
slope  can  be  summed,  since  the  analysis  is  linear. 

This  method  applies  to  flow  over  terrain  which  does  not  possess  steep  slopes,  for 
meteorological  conditions  that  can  be  characterized  by  relatively  simple  profiles 
of  windspeed  and  temperature  in  the  vertical.  If  the  stratification  is  strong 
(that  is,  the  Froude  number  u/NH  <  1)  these  methods  are  not  applicable.  As 
discussed  in  connection  with  the  diagnostic  models,  the  flow  tends  to  be  two- 
dimensional  in  a  layer  below  the  critical  di.viding- streamline .  This  means  that 
perturbations  are  as  large  as  the  incident  flow  in  this  layer,  which  invalidates 
the  model.  For  general  application,  the  FLOWSTAR  approach  would  have  to  incorpo¬ 
rate  a  complementary  model  in  order  to  treat  this  aspect  of  the  flow.  Note, 
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however,  that  the  flow  above  the  critical  dividing-streamline  height  behaves  as 
a  weakly  stratified  flow,  so  that  the  outer  region  solutions  might  be  adapted  to 
model  this  region. 


Figure  2.  Flow  regions  for  the  linear  analysis  (Carruthers  and  Hunt,  1990). 


2.1.4  Prognostic  flow  models 

Prognostic  models  used  to  simulate  mesoscale  flows  in  complex  terrain  were 
recently  reviewed  by  Pielke  (1989),  and  summarized  by  Pielke  et  al.  (1990). 
Seven  distinct  models  developed  and  applied  in  North  America  were  included  in  the 
review.  Of  these,  four  were  found  to  be  applicable  to  complex  terrain- -the  MASS 
model,  the  NCAR/Penn  State  (MM4)  model,  the  CSU  RAMS  model,  and  the  Los  Alamos 
model  HOTMAC  (Yamada  and  Bunker,  1988).  Each  of  these  includes  a  primitive 
equation  framework,  a  nested  grid,  and  each  makes  use  of  the  hydrostatic 
approximation.  Only  CSU  RAMS  also  includes  a  nonhydrostatic  option.  The  minimxim 
spacing  allowed  for  the  horizontal  grid  varies  among  these  models.  CSU  RAMS  can 
accept  a  spacing  of  about  100  m;  the  Los  Alamos  model  has  a  minimum  spacing  of 
380  m;  MM4.  has  a  minimum  of  about  1000  m;  and  MASS  requires  a  horizontal  spacing 
of  at  least  7000  m. 

The  four  models  that  are  readily  applicable  to  complex  terrain  include,  in 
principle,  a  complete  description  of  the  mean  flow  field,  including  diabatic 
effects.  However,  by  their  very  nature,  they  require  significant  computer 
resources.  Most  simulations  make  use  of  super  computers,  but  "super-mini"  or 
graphics  super-workstations  appear  to  be  viable  alternatives.  In  a  recent 
demonstration  project.  Mattocks  et  al.  (1989)  adapted  MASS  for  execution  on  an 
Alliant  FX-1  minisupercoraputer ,  with  satisfactory  results.  Their  tests  indicate 
that  the  minisupercoraputer  ran  at  approximately  l/20th  the  speed  of  a  CRAY  2 
supercomputer,  which  is  about  three  times  the  anticipated  speed. 

2.2  Turbulence  Characterization  Approaches  and  Parameterizations 

The  various  flow  field  modeling  approaches  outlined  in  section  2.1  characterize 
turbulence  in  the  flow  in  several  ways.  The  objective  analysis  and  the  diagnos¬ 
tic  methods  do  not  rely  on  characterizing  the  turbulence  at  all,  and  at  most 
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include  estimates  of  turbulence  in  an  ad-hoc  manner  if  diffusion  estimates  are 
required.  In  such  a  case,  the  turbulence  may  simply  be  parameterized  in  terms 
of  a  Pasquill  stability  class  index,  or  it  may  be  explicitly  related  to  micro- 
meteorological  parameters.  The  point  is  that  the  turbulence  is  derived  from  mean 
properties  of  the  flow  field- -the  flow  field  itself  is  determined  without  regard 
for  the  effect  of  turbulence. 

Analytic  flow  field  models  do  incorporate  aspects  of  the  turbulence  in  modeling 
the  mean  flow.  This  is  done  primarily  in  the  inner  region,  where  turbulence  is 
in  local  equilibrium.  Essentially,  both  the  shear  stress  due  to  the  mean  flow 
and  the  perturbation  shear  stress  are  explicitly  contained  in  the  linearized 
equations  of  motion.  These  shear  stresses  are  typically  characterized  in  terms 
of  a  mixing  length  formulation. 

The  prognostic  models  also  explicitly  incorporate  models  for  turbulent  transport, 
which  make  up  an  integral  part  of  the  system  of  equations  used  to  obtain  the  mean 
flow  field.  Depending  on  the  type  of  closure  employed,  the  turbulence  is  modeled 
in  terms  of  a  gradient  transfer  with  specified  eddy  coefficients  (K-Theory) ,  or 
by  formulating  the  Reynolds  stresses  explicitly,  and  closing  the  system  of 
equations  with  a  model  for  the  higher-order  products  of  the  turbulent  velocity 
fluctuations  (for  example,  second-order  closure).  The  MASS  and  NCAR/Penn  State 
models  use  first-order  closure,  the  Los  Alamos  model  uses  second-order  closure, 
and  the  CSU  RAMS  model  includes  an  option  to  use  either . 

Other  efforts  to  model  atmospheric  turbulence  are  found  in  the  literature, 
although  these  have  not  been  coupled  to  wind  field  models  suitable  for  applica¬ 
tions  in  complex  terrain.  Two  such  approaches  are  LES  (Large  Eddy  Simulation) 
and  KS  (Kinematic  Simulation).  LES  is  an  approach  in  which  a  modified  form  of 
the  Navier-Stokes  equation  is  solved  on  a  grid  sufficiently  small  to  resolve 
the  eddies  of  interest.  Subgrld-scale  motions  are  parameterized,  and  must 
generally  be  small  enough  to  lie  with  the  inertial  subrange.  LES  models  similar 
to  Deardorff  (1974)  for  the  convective  boundary  layer  use  gradient  transfer 
equations  for  subgrid-scale  exchange,  and  include  an  equation  for  the  energy  of 
these  motions.  As  a  result,  LES  results  for  the  stable  boundary  layer  may  be 
significantly  influenced  by  the  particular  parameterization  used  for  subgrid 
scale  processes. 

KS  is  a  recently  proposed  method  for  simulating  the  details  of  a  turbulent 
velocity  field  with  prescribed  spectral  properties.  As  KS  theory  begins  with  an 
analytic  Fourier  representation,  it  is  not  tied  to  a  computational  grid. 
However,  it  contains  no  dynamics,  so  that  the  effects  of  vortex-  stretching,  for 
example,  are  represented  only  to  the  extent  that  the  spectral  properties  embody 
the  result  of  this  process.  Because  of  this,  the  KS  approach  may  be  viewed  as 
a  "repository"  of  present  knowledge  about  Eulerian/Lagrangian  statistics  of 
turbulence.  KS  theory  will  be  examined  in  greater  detail  in  section  2.4. 

In  the  remainder  of  this  section,  we  present  methods  for  characterizing  turbu¬ 
lence  variances  (as  and  for  example)  _n  undisturbed  flows  by  means  of 
established  micrometeorological  parameters,  and  we  present  related  expressions 
for  time  scales.  Approximate  expressions  for  the  spectra  of  turbulence  are  also 
discussed.  Because  flows  are  altered  by  the  presence  of  terrain,  properties  of 
the  turbulence  are  altered  as  well,  and  suggestions  for  how  to  formulate  these 
changes  are  reviewed. 
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2.2.1  and  a. 

Many  laboratory  experiments,  field  studies,  and  numerical  simulations  (for 
example,  Deardorff  and  Willis,  1975;  Caughey,  1981;  Lamb,  1981)  have  shown  the 
importance  and  utility  of  convective  scaling  in  the  convective  boundary  layer. 
Convective  scaling  has  been  successfully  applied  to  data  collected  at  a  wide 
variety  of  sites,  including  oceans,  rural  land  surfaces  (for  example,  Hicks, 
1985)  and  urban  areas  (Ching  et  al.,  1983).  Similarly,  in  the  stable  boundairy 
layer,  local  scaling  has  been  shown  to  apply  (for  example.  Hunt,  1981; 
Nieuwstadt,  1984). 

Weil  (1985)  and  Briggs  (1985)  provide  reviews  on  the  use  of  similarity  theory  in 
diffusion  models.  In  the  convective  boundary  layer,  Weil  describes  the  turbu¬ 
lence  characteristics  in  three  layers : 

(1)  Surface  layer  -  z  <  0.1  h;  ~  constant  with  height 

increases  with  height 

(2)  Mixed  layer  -  O.lh  <  z  <  0.8h;  -  constant  with  height 

-  constant  with  height 

(3)  Entrainment  layer  -  z  >  0.8h;  decreases  with  height 

decreases  with  height 

In  the  surface  layer,  Panofsky  et  al.  (1977)  propose  the  following  relations. 

«  U,  [4  0 .6 


a^=u,[l.6  +  2.9 


where  u*  is  the  surface  friction  velocity  (m/s) ,  and  L  is  the  Monin-Obukhov 
length  (m) . 

Hicks  (1985)  suggests  the  following  for  the  mixed  layer  (0.1  to  0.8  h) . 

=  (3 .6  U.2  +  0.35 


=  (1.2  U.2  +  0.35  C5) 
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In  the  neutral  boundary  layer,  Arya  (1984)  reports  monotonically  decreasing 
values  of  cr„  and  throughout  the  mixed  layer.  Using  the  Blackadar  and  Tennekes 
(1968)  relationship  for  the  neutral  boundary  layer  height,  Arya's  results  can  be 
expressed  as; 

<Jy  =  1.5  exp  {-0.9  z/h)  (6) 

=  1.3  exp (-0.9  z/h)  O) 


In  the  stable  boundary  layer,  Nieuwstadt  (1984)  finds  that  and  cr„  bear 
constant  ratios  with  the  local  friction  velocity. 

Oy/u.,  =  Cy  (8) 

(9) 


where  u*,  is  the  local  friction  velocity  (m/s),  and,  Cv  and  are  constants. 

Hanna  et  al.  (1986)  suggest  that  C,  »  1.6.  C.  has  a  value  =«  1.3  (Nieuwstadt, 
1984).  The  local  friction  velocity,  u*, ,  can  be  expressed  (Nieuwstadt,  1984)  as: 

U.,  =  U.  (1  -  (10) 


2.2.2  Lagrangian  time  scale 

Empirical  equations  for  the  Lagrangian  time  scale  have  been  proposed  by  Hanna 
(Nieuwstadt  and  Van  Dop,  1982)  for  use  in  Monte  Carlo  particle  modeling.  He  also 
presents  expressions  for  <Tv  and  <t„,  and  for  a^.  Although  these  expressions 
differ  from  those  presented  above,  we  point  out  that  he  equates  cr„  with  for 
the  convectively  unstable  boundary  layer,  and  equates  cr„  with  for  the  neutral 
and  stable  stratified  boundary  layer.  For  the  latter,  is  given  as: 

=  2.0  u*  /l  *  -^1  {stable)  (11) 


dy  =  2.0  u.  exp (-3  t  z/u,)  {neutral)  (12) 
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where  f  is  the  coriolis  parameter.  The  Lagrangian  time  scales  proposed  by  Hanna 
are: 


Convectivelv  Unstable  Boundary  Laver 

=  0,15  h/a^ 


(13) 


Z  <  O.lA;  z  -  Zg,  >  -L: 


0.1  z/a^ 

0.55  +  0.38  (z  -  Zj/L 


(14) 


<  O.IA;  z  -O  <  -L;  =  0.59  z/o^ 


(15) 


z  >  O.lA:  =  0.15  |l  -exp  |-5-| 


(16) 


Stably  Stratified  Boundary  Laver 


=  0.15 


z|0-s 

OuU/ 


(17) 


h  /  ^  \  0  •  5 

2V  =0.07—4 
o^\hl 


(18) 


=0.10 


I  Z\®'® 

(a) 


(19) 


Neutral  Boundary  Laver 


=■  = 


0.5  z/o^ 

1  +  15  f  z/u. 


(20) 


In  these  expressions,  the  subscript  "a"  denotes  a  quantity  describing 
fluctuations  along  the  flow,  and  "c"  denotes  a  quantity  describing  crosswind 
fluctuations.  We  use  the  alternate  notation  in  which  "u"  is  along  the  flow,  and 
"v"  is  across  the  flow  (horizontal)  in  describing  ay  and  and  for  o-„. 
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2.2.3  Velocity  spectra 

Following  the  work  performed  by  Kaimal  et  al.  (1972,  1976)  and  Kaimal  (1973)  in 
characterizing  velocity  spectra,  several  authors  have  developed  revisions  which 
extend  the  expressions  to  greater  heights  above  the  surface.  For  the  convective 
boundary  layer,  Hojstrup  (1982)  has  developed  expressions  that  are  valid 
throughout  the  lower  half  of  the  layer,  and  that  smoothly  approach  the  neutral 
behavior  found  in  the  Kansas  data  (Kaimal  et  al.,  1972).  More  recently,  Moraes 
(1988)  extended  Kaimal' s  (1973)  formulas  for  the  stably  stratified  surface  layer 
to  the  upper  part  of  the  boundary  layer  by  incorporating  Nieuwstadt's  (1984) 
local  similarity  scaling.  The  remaining  area  in  which  the  representations  may 
be  inadequate  is  the  upper  half  of  the  convective  boundary  layer.  Hojstrup 
(1982)  points  out  that  spectra  in  this  region  are  influenced  by  entrainment 
processes  at  the  top  of  the  layer,  and  the  effect  of  the  "lid"  itself  as  a 
boundary  that  limits  the  vertical  scale  of  motion. 

2. 2. 3.1  Convective  boundary  layer  (Hojstrup,  1982) 


nS^(n) 


1  +  2.2  \  -L/ 


(1  +  33  (1  +  15z/A) 


nS^(n) 


(1  +  2/i)5/3  \  -LJ 


_ 17  4^(1  -  z/A)2 _ 

(1  +9. 54^)  5/3(1  +  2.8z/73)2/3 


nS^(n) 

=  F{f,z/h) 


0.95fi 


(1  +  2fi)5/3 


4\2/3  ^  2f(l  -  Z/A)2 
1+5.3  (f)5/3 


(23) 


where  n  is  the  frequency  in  Hz  and  the  reduced  frequencies,  f,  are  given  as 


f  =  nz/u  ;  -  f/(l  +  15z/h) 

=  nh/u  ;  fj.^  =  f/{l  +  2.8z/A) 


f{f,  z/h) 


f2  >  (0.3z/A)2 

+  (0.15)2 


(24) 
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Also,  h  is  the  height  of  the  lowest  inversion,  L  is  the  Monin-Obukhov  length,  and 
is  the  friction  velocity  in  the  surface  layer. 

Such  analyses  based  on  observed  data  from  the  convective  boundary  layer  can  be 
significantly  augmented  by  simulation  data  bases  produced  by  LES  models.  The 
data  base  and  analyses  performed  by  Moeng  and  Wyngaard  (1988)  should  greatly 
facilitate  a  comprehensive  modeling  of  the  CBL. 

2. 2. 3. 2  Neutral/stable  boundary  layer  (Moraes,  1988) 

The  velocity  spectra  in  the  surface  layer  follow  Kaimal  (1973) ,  as  recast  by 
Moraes  and  Epstein  (1987) : 

=  V  0.164/  ,25) 

^  0.164/5/3 


where  i  denotes  one  component  of  the  turbulence  (u,v,or  w) ,  f  -  nz/u,  and  f  is 
the  dimensionless  dissipation  rate  for  turbulent  energy.  The  parameter  has 
the  following  values : 

Tu  -  0.3 

7v  “  0.4 

-  0.4 

The  reference  value  for  the  frequency  scale  is  given  by 

^Oi  =  Pifl  ^  2.5(z/i)5/5]3/2  (26) 

where  i  again  denotes  one  of  the  components  u,  v,  or  w.  The  parameter  has  the 
following  values : 

-  0.012 
0^  -  0.045 
0^  -  0.094 

Above  the  surface  layer,  Moraes  (1988)  retains  the  functional  form  of  spectra  in 
the  surface  layer,  but  incorporates  results  from  local  similarity  theory.  The 
resulting  equation  for  the  velocity  spectra  is 

nSjin)  ^  _ 0.164/ _  ,27) 

+  0.164/5/3 
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where  U,  is  the  local  turbulence  velocity  scale,  is  the  reference  value  for 
the  frequency  scale  in  the  stable  boundary  layer,  is  the  dimensionless  energy 
dissipate  rate  in  the  stable  boundary  layer: 

»  1.2(1  ♦  3.7Z/A)  (28) 


and  7i  and  f  are  the  same  as  in  the  surface  layer.  The  local  scaling  length  A 
is  related  to  the  stable  boundary  layer  (h)  and  the  Monin-Obukhov  length  (L)  for 
the  surface  layer.  Following  Nieuwstadt  (1984),  Horaes  finds  that 


A  =  (1  -  L 


(29) 


Furthermore,  the  frequency  scale  for  the  stable  boundary  layer  is  given  by 


(30) 


where 


Q  5B£.  -  A  [l  -H  0.628 

1  +  0.422  (h/L) 


(31) 


Values  for  the  parameter  are  the  same  as  those  in  the  surface  layer. 

The  recent  data  reported  by  Smedman  (1991)  should  also  improve  characterization 
of  the  spectra  in  the  SBL. 

2.2.4  Changes  Induced  by  terrain 

Carruthers  and  Hunt  (1990)  provide  an  overview  describing  how  changes  in  the  mean 
flow  over  hills  also  affects  the  properties  of  turbulence.  In  essence,  they 
state  that  the  turbulence  is  in  near  equilibrium  with  the  altered  flow  in  the 
inner  region  near  the  surface  of  the  hill,  and  the  turbulence  in  the  outer  region 
is  changed  by  the  distortion  of  the  mean  flow.  The  two  regions,  inner  and  outer, 
largely  conform  to  the  inner  and  outer  regions  referred  to  in  figure  2. 

In  the  inner  region,  the  turnover  time  scale  for  the  eddies  is  of  order  z/u. , 
which  is  less  than  the  travel  time  over  the  hill.  Consequently,  the  turbulence 
in  the  flow  has  time  to  adjust  to  changes  in  the  flow,  and  is  therefore  in 
approximate  equilibrium  with  the  flow.  Therefore,  fractional  changes  in  the 
turbulence  are  equal  to  fractional  changes  in  the  shear  stress: 


^  =  Al  (32) 
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Because  the  change  in  the  shear  stress  decreases  rapidly  with  height  above  the 

surface,  it  is  found  that  the  Auf  decreases  t  nearly  zero  at  one-third  the 
depth  of  the  inner  region. 

The  turnover  time  increases  with  the  size  of  the  eddies  and,  therefore,  with 
height,  so  that  it  exceeds  the  travel  time  over  the  hill  for  eddies  in  the  outer 
layer.  Here,  changes  in  the  flow  field  distort  the  eddies  in  the  approach  flow, 
and  the  effect  can  be  modeled  by  rapid  distortion  theory  as  explained  by  Britter 
et  al,  (1981).  Carruthers  and  Hunt  (1990)  report  that  Newley  (1985)  and  Zeman 
and  Jensen  (1987)  have  found  that  anistropy  in  the  upstream  turbulence  should  be 
included,  and  suggested  that 


)  AS 


(33) 


4^  =  .6(1  -  AS  (34) 

uj  6  3 


where  AS  is  the  speed-up  factor.  This  result  applies  for  the  special  case  in 

■■X 

which  uj  =  uj,  where  component  1  is  along  the  flow,  2  is  across  the  flow 
(lateral),  and  3  is  across  the  flow  (normal  to  the  terrain).  Finnigan  (1988) 
elaborates  on  the  processes  that  are  important  between  the  inner  and  outer 
regions,  which  include  curvature  effects  and  nonlinear  interactions  as  well  as 
rapid  distortion. 

2.3  Overview  of  Monte  Carlo  Lagrangian  Particle  Modeling 

In  Lagrangian  particle  models ,  particles  are  moved  in  the  Eulerian  frame  by  the 
basic  equation 


2:(t  +  At)  =  jc(t)  -t-  t)  +  jZ  (j:,  t)  ]  •  A  t 


(35) 


where  u(x,t)  represents  the  space-time  varying  mean  flow  velocity,  and  u'(x,t) 
is  the  stochastic  vector  component.  Evolution  of  this  stochastic  component  is 
then  modeled  as  a  first-order  autoregressive  or  Markov  process  using  Langevin 
stochastic  differential  equations  (Gifford,  1982;  Sawford,  1984), 

du'  =  -u'dtlTj^  +  du" ,  (36) 
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for  each  vector  component  and  where  is  the  Lagrangian  time  scale  and  du"  are 
the  random  velocity  components  arising  from  a  white  noise  modeling  of  particle 
accelerations.  Sawford  (1991)  discusses  equation  (36)  in  relation  to  lower  and 
higher  order  autoregressive  equations  involving  the  particle's  position  (chat  is. 
random  walk  modeling)  and  acceleration,  respectively. 

The  stochastic  term  u"  is  often  reexpressed  (for  example,  Rodean,  1991)  as 
(C(je)l/2dW(t)  ,  where  is  the  universal  constant  for  the  Lagrangian  structure, 
e  is  the  turbulent  kinetic  energy  dissipation  rate  and  dW  is  an  increment  of  a 
Wiener  process  having  zero  mean  and  a  variance  of  dt. 

In  many  applied  models,  the  solution  of  equation  (36)  is  expressed  as 


u'(t)  *  ■  uUt  -  At)  ♦  i?  •  Ou  •  y/l  -  -f  du 


(37) 


where  -  exp(-At/TLu) 

Tu  is  the  Lagrangian  time  scale  for  the  u  component, 

R  is  a  random  Gaussian  number  having  zero  mean  and  unit  variance, 

<7u  is  the  standard  deviation  of  u  velocity  fluctuations,  and 
dy  is  the  u  component  drift  velocity. 

Similar  equations  are  used  for  the  v'  and  w'  wind  components  as  well.  The  drift 
velocity  can  be  shown  to  be  a  necessary  ingredient  to  prevent  buildup  of 
particles  in  low  turbulence  regions  and  is  specified  by  Sawford  (1985)  as 


du  =  (1  -  ttu)  • 


dx 


[  (0^2  +  u2)/(20u2)  ] 


(38) 


for  the  u  component.  Other  authors  (for  example,  Legg  and  Raupach,  1982;  Ley  and 
Thompson,  1983)  model  this  drift  velocity  as  effectively,  but  more  simply,  as 


(39) 


Equation  (37)  reduces  to  the  random  walk  model  in  the  limi  and  Co  a 
purely  deterministic  model  as  Tl^  -►  «.  Further,  a  moments  analysis  of  this 
equation  (37)  shows  that  while  it  has  a  mean  corresponding  to  the  drift  velocity, 

dy,  and  variance  corresponding  to  (l  -  the  third  moment  (or  skewness) 
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vanishes.  Thus,  equation  (37)  is  not  directly  appropriate  to  the  modeling  of 
convective  conditions. 

The  somewhat  more  complex  problem  of  skewed  turbulent  spectra  is  basically  solved 
by  selecting  the  random  component,  u" ,  of  equation  (36)  from  the  desired 
probability  density  function  (pdf)  of  turbulent  velocities .  Weil  (1990)  provides 
a  comprehensive  review  of  such  approaches  as  applied  to  modeling  the  convective 
boundary  layer. 

Another  issue  requiring  modification  of  the  Langevin  formulation  involves 

attempts  to  generate  finite  velocity  cross  correlations  (for  example,  <  0) 

rather  than  the  zero  cross  correlations  which  result  from  the  evolution  of 
independent  Langevin  equations  for  the  u' ,  v' ,  and  w'  components.  Zannetti 
(1990)  provides  a  review  of  such  attempts  and  has  included  such  cross 
correlations  into  his  MC-LAGPAR  II  model  (Zannetti,  1986).  In  our  planned  usage 
of  Langevin  equation  models  for  only  the  shorter  wavelength  portion  of  the 
turbulent  spectrum,  it  will  not  be  necessary  to  add  such  cross  correlation 
contributions  as  they  will  arise  naturally  from  the  flow  fields  specified  via 
kinematic  simulation. 

Finally,  we  note  that  equation  (35)  is  strictly  appropriate  for  massless  and 
inertialess  particles  that  instantly  adjust  to  the  turbulent  velocities.  For 
particles  with  significant  gravitational  settling  velocities,  Vg,  or  inertia, 
equation  (35)  is  replaced  by  the  pair  of  vector  equations 

jC<t  +  At)  =x(C)  *  y(jc,  t)  At  (^0) 

and 

=  yiuix,  t)  *  uHx,  t)  -  c))  -  s 


where 

v(2c,t)  is  the  velocity  of  the  particle, 

g  -  (0,0, g)  where  g  is  the  acceleration  due  to  gravity,  and 

7  “  1/^p  is  the  inverse  of  the  particle  response  time  Tp. 

Clearly  as  the  particle  response  time  goes  to  zero  equation  (41)  has  the 
solution,  Y  —  u  +  u' ,  and  equation  (35)  is  recovered.  For  spherical  particles 
governed  by  Stokes  law.  Squires  and  Eaton  (1991)  assume 

7  »  IQp/ (Ppd^)  (^2) 
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where  n  is  the  dyiianiic  viscosity  of  air,  Pp  is  the  particle  density,  and  d  is 
particle  diameter.  Mallier  and  Maxey  (1991)  discuss  the  much  more  general 
equations  applicable  to  nonspherical  particles.  Frost  et  al .  (1982)  use  a 
potentially  more  generally  expression  for  their  MoCaPD  model  which  reduces  to 

Y  =  I  (p^/pp^Cd 


where  is  the  density  of  air, 

is  the  drag  coefficient  defined  as  a  function  of  particle  Reynolds 
number ,  and 

|Vj|  is  the  magnitude  of  the  relative  velocity  of  the  particle  in  air. 

However,  a  recent  study  by  Hashem  and  Parkin  (1991)  suggests  that  the  effects  of 
particle  inertia  are  negligible  for  particles  up  to  500  /im  in  diameter  so  that 
1  can  be  expressed  in  terns  of  the  terminal  velocity,  v,.,  of  the  particle  and  the 
turbulence  parameters  and  Tl„  for  the  w  component.  This  then  leads  to  a 
closed  form  expression  for  the  particle  velocity. 

2.4  Kinematic  Simulation  and  Hybrid  Theories 

In  the  case  of  single  particle  Lagrangian  trajectories,  individual  particle 
velocity  components  (u^',  Vj_',  Wj^')  evolve  independently  of  one  another  and  are 
uncorrelated  with  the  velocity  components  of  a  neighboring  particle.  These 
characteristics  of  the  particles'  motions  imply  that  the  flow  fields  influencing 
a  specific  particle  are  not  assured  of  being  divergence  free  and  that  colocated 
particles  are  moving  completely  independently  of  one  another.  The  first  of  these 
deficiencies  can  lead  to  the  accumulation  of  particles  in  low  turbulence  zones: 
a  problem  that  can  be  compensated  for  on  average  by  introducing  drift  velocities. 
The  second  deficiency  implies  that  such  single  particle,  Monte  Carlo  Lagrangian 
trajectory  modeling  can  only  lead  to  ensemble  average  predictions  and  statistics. 

The  most  logical  way  to  avoid  these  compromises  is  to  let  the  particles  follow 
the  actual  turbulent  flow  fields  present  at  any  point  in. space  and  time.  Such 
"actual"  turbulent  flow  fields  can  either  be  generated  by  dynamical  simulations 
(DS) ,  involving  solutions  of  the  Navier  Stokes  equations  over  space  and  time,  or 
by  kinematic  simulations  (KS) ,  involving  the  conjecture  of  flow  fields  which  at 
a  minimiim  obey  the  constraint  V-V  -  0  (or  -  0). 

DS  methods  encompass  both  the  direct  numerical  simulation  (DNS)  and  large  eddy 
simulation  (LES)  methodologies  and  lead  to  a  more  complete  description  of  the 
turbulent  flow  fields  in  the  sense  that  all  (or  most)  of  the  relevant  physics  is 
included.  Such  approaches  are  currently  too  computationally  intensive  for  the 
intended  application  but  can  be  utilized  as  an  important  resource  to  guide  in  the 
development  of  appropriate  KS  field  properties  and  supplemental  constraints. 

KS  theory  begins  with  the  representation  of  the  vector  flow  field  u(x,t)  as  a 
four- dimensional  Fourier  transform  of  an  arbitrary  amplitude  function,  S(k,u)), 
where  k  is  the  wavenumber  vector  and  w  represents  frequency.  This  continuou.s 
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transform  is  then  discretized  into  finite  sums  following  the  work  of  Kraichnan 
(1970)  and  Drummond  et  al.  (1984).  Further  assuming  that  a  relatively  narrow 
range  of  w  contributes  at  a  particular  k  for  the  larger  eddies,  Fung  et  al . 
(1490)  develop  the  representation  of  the  sub  -  grid- scale  flow  field  u’(x,c)  in 
terms  of  N  spectral  components  as 


u' (j£,  t)  =  ]cosag(t) 

n  (44) 

+  1^(  t)  ®  ^  ]  3in  ap(  t)  } 
n 


where 


ao{  t) 


-  lit 


t 


(45) 


u  is  the  mean  flow, 


3  -  At)  +  u„„(t) 


(46) 


At  represents  a  migrating  spatial  offset  position  of  the  nth  spectral 
component, 

and  the  vector  amplitudes  a„{ t)  ,  b„{  t) 

and  eddy  migration  velocity  ^  ( t) 

—  n 

each  evolve  according  to  a  Langevin  type  equation  with  a  relevant  Lagrangian  time 
scale , 

T,  =  l/a>„  =  i/(a„.^  -k^),  (47) 


with  representing  the  velocity  variance  of  the  nth  spectral  component.  As 
the  components  of  u’(x,c)  are  generated  by  vector  cross  products  (that  is,  ®)  of 

the  a,  b  amplitudes  and  Che  ic^j  unit  vector,  the  divergence  relation,  Vou'  -  0, 
is  equivalent  to  the  sura  on  n  of  the  dot  products,  is  therefore 
immediately  satisfied. 

Equation  (44)  therefore  contains  nine  randomly  evolving  quantities  for  each 
spectral  component  n.  The  fact  chat  the  sine  and  cosine  pieces  evolve 
independently  implies  that  neither  standing  wave  nor  traveling  wave  (chat  is, 
6xp[±  iOo(t)])  solutions  are  allowed  to  persist.  The  definition  of  the 
appropriate  time  scale,  T^,  and  the  practical  constraint  chat  Langevin  equations 
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should  be  marched  at  a  time  step.  At,  such  that  At  <*  T„,  implies  that  the 
smallest  eddies  modeled  (that  is,  largest  k„)  will  dictate  the  time  step  for 
advancing  particle  trajectories.  Among  other  problems,  excessively  large  time 
steps  would  imply  that  particles  would  not  follow  local  streamlines  and  a 
principal  feature  of  the  divergence  free  flow  aspect  of  the  formulation  (that  is, 
avoidance  of  particle  accumulation  in  zones  of  low  turbulence)  would  be  lost. 
This  constraint  on  At  suggests  that  the  turbulent  spectrum  modeled  by  the  KS 
approach  be  chopped  off  at  a  k  value  representing  an  appropriate  trade  off  in 
realism  and  computational  cost,  and  the  effect  of  the  high  frequency  remainder 
of  the  spectrum  be  included  as  a  random  walk  or  Langevin  equation  component  of 
the  particle's  motion.  Such  an  approach  has  recently  been  suggested  in  Sakai  et 
al.  (1991),  by  the  Cambridge  University  team,  responsible  for  many  of  the 
advances  in  KS  theory,  as  a  more  economical  alternative  to  their  original 
approach  (for  example,  Fung  et  al.,  1990)  involving  separate  but  explicit 
treatment  of  the  smaller  "inertial"  eddies  that  are  carried  along  by  the  larger 
"sweeping"  eddies  considered  here. 

Fung  et  al.  (1990)  also  provide  guidance  on  developing  vector  amplitudes, 
(a„,  b„)  ,  the  random  component  of  the  eddy  migration  velocity,  U„„,  and  choice 
of  appropriate  wave  number  values  Iq,.  Each  of  the  above  9N  components  of  a,  b, 

and  are  chosen  to  have  zero  mean,  a  variance  of  oV a'  absence  of 
covariance  terms.  Despite  this  lack  of  covariance  in  the  vector  amplitudes, 
turbulent  velocity  components  will  display  appropriate  negative  normalized 
covariances  of  order  -  1/2, 

The  equation  (44)  representation  of  the  turbulent  flow  field  is  most  appropriate 
for  isotropic  turbulence.  In  the  case  of  anisotropic  turbulence,  the  time 
scales,  T„,  of  equation  (47)  could  easily  be  made  vector  quantities;  however,  due 
to  the  generation  of  velocity  components  via  a  cross  product  relation  (for 
example,  the  x  component  of  u'  involves  both  the  y  and  z  components  of  the  a,  b 
amplitudes),  it  becomes  difficult  to  drive  the  independent  Langevin  equations  for 
the  amplitudes  with  an  appropriate  time  scale.  In  addition,  equation  (44)  is 
most  appropriate  for  an  unbounded  medium  (that  is,  no  obstacles).  The  presence 
of  a  solid  surface,  such  as  the  ground  at  z  -  0  may  be  handled  by  the  method  of 
images  or  may  be  designed  explicitly  for  that  geometry,  as  will  be  done 
subsequently  for  a  surface  at  z  -  0,  In  the  case  of  the  surface  at  z  -  0,  the 
image  method  would  suggest  the  combination 

odd  ^  [u' Z,  t)  -  U' {X,y,  ~Z.  t)] /2  (48) 

as  that  appropriate  for  ensuring  that  only  those  velocity  components  that  are  odd 
in  z  would  suir/ive .  This  procedure  will,  of  course,  kill  off  all  threee 
components  of  the  resultant  u' odd  at  z  -  0,  whereas  it  may  be  desired  to 
eliminate  only  the  z  component  of  velocity  at  the  surface. 

The  limitations  of  turbulence  Isotropy  and  suppression  of  the  normal  component 
(w')  at  a  z  -  0  surface  that  are  embedded  within  equation  (48)  may  be  eliminated 
by  developing  a  customized  set  of  equations  for  the  turbulent  velocity 
components.  One  such  possibility,  for  the  nth  spectral  component,  is  given  as 


v/  =  [33  cos  +  233  sin  a3]  sinCic^  •  z) 

u'  =  (k,/kx)  [-ajP^  sin  a,  +  Jbj  y^cos  03] co,  ‘ 

+  {ky/k(.)  [a^  sin  cos  cos  sin  Oj] 

v'  =  ikjky)  [-a3(l-p2)  sin  a,  +  b-^  (1-Y^)  cos  a3]  cos  {kg  •  z) 

+  {k^/k^)  [-ay  cos  sin  -  by  sin  cos  ttj] 

where 

Qti  -  ky(y  -  Yoi  -  vt)  +  k,(z  -  Zqi  -  wt)  , 

02  -  k,(x  -  Xo2  -  ut)  +  k,(z  -  Zo2  *  wt)  , 

ag  -  k,(x  -  Xo3  -  ut)  +  ky(y  -  y^g  -  vt)  , 

kj  -  (k,^  +  ky^)',  and 

the  four  a,  b  amplitudes,  the  six  coordinates  (Xo2.  x^g,  yoi,  yog,  Zqi,  Zqz)  the 
two  projection  operators  P,  7  are  all  allowed  to  evolve  in  time  via  separate 
Langevin  equations.  It  should  also  be  noted  that  the  nearly  ubiquitous  subscript 
n  has  been  dropped. 

In  the  above  equations,  the  (ag,  bg)  amplitudes  determine  the  strengths  of  two- 
dimensional  eddies  lying  in  the  x-z  and  y-z  planes  with  the  projection  operators 
(^.7)  defined  such  that  -1  <  /3,  7  <  +1,  determine  the  apportionment  between  x  and 
y  fractions.  Similarly  the  (a^,  b^)  amplitudes  determine  the  strengths  of  two- 
dimensional  eddies  and  convergence/divergence  zones  in  the  x-y  plane. 

A  principal  advantage  of  the  equation  (49)  formulation  is  that  the  time  scales, 
Tg,  which  govern  the  evolution  of  the  vertically  oriented  eddies  can  be 
completely  decoupled  from  the  time  scales,  Ty,  of  the  longer  lived  transverse 
fluctuations  in  the  x-y  plane  given  by  equations  (50)  and  (51).  We  also  note 
that  the  simple  z  dependence,  sin(kj«z),  could  easily  be  generalized  to  a  more 
general  form,  f(z) ,  with  the  additional  requirement  being  that  the  cos(kj*z) 
dependence  in  u'  ,  v'  be  replaced  with  df(z)/dz.  Additional  z  dependence  may  also 
be  added  to  the  transverse  eddy  components  as  needed  to  provide  more  realistic 
vertical  profiles  of  turbulent  velocity  moments. 

3.  MULTI -SPECTRAL  SMOKE  ASPECTS  OF  THE  PROBLEM 

3.1  Dry  Deposition  of  Particles 

Many  complex  processes  are  involved  in  the  transfer  and  deposition  of  pollutants 
at  the  surface.  Sehmel  (1980)  compiled  a  list  of  the  factors  known  to  influence 
dry  deposition  rates  (see  table  1).  The  variables  listed  include  the  properties 
of  the  depositing  material  (for  example,  particle  size  and  density;  gas 
diffusivity,  solubility,  and  reactivity),  the  characteristics  of  the  surface  (for 
example,  surface  roughness,  vegetation  type,  amount,  and  physiological  state), 
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and  atmospheric  variables  (for  example,  stahility,  turbulence  intensity).  Hicks 
(1982)  noted  the  important  differences  controlling  the  deposition  of  larger 
particles  (for  example,  gravitational  settling,  inertial  impaction)  and  those 
controlling  gases  (for  example,  turbulence,  molecular  diffusion).  Deposition  of 
small  particles  is  complicated  by  the  fact  that  they  may  be  influenced  by  the 
processes  affecting  both  gases  and  large  particles.  Figure  3  shows  the 
tremendous  variability  of  deposition  velocity  to  particle  size  alone. 

TABLE  1.  FACTORS  INFLUENCING  DRY  DEPOSITION  RATES* 


Micrometeorological 
_ Variables _ 


Depositing 

Material 


Surface 

Variables 


ParticLes. 


Aerodynamic  roughness 

-  Mass  transfer 

(a)  Particles 

(b)  Gases 

-  Heat 

•  Momentum 

Atmospheric  stability 
Diffusion,  effect  of; 

-  Canopy 

-  Diurnal  variation 

-  Fetch 

Flow  separation: 

-  Above  canopy 
'  Below  canopy 

Friction  velocity 
Inversion  layer 
Pollutant  concentration 
Relative  humidity 
Seasonal  variation 
Solar  radiation 
Surface  heating 
Temperature 
Terrain 

-  Uniform  ■ 

-  Nonuniform 
Turbulence 
Wind  velocity 

Zero -plane  displacements 

-  Mass  transfer 

(a)  Particles 

(b)  Gases 

-  Heat 

-  Momentum 


Agglomeration 

Diameter 

Density 

Diffusion 

-  Brownian 

-  Eddy  equal  to 

(a)  Particle 

(b)  Momentum 

(c)  Heat 

'  Effect  of  canopy  on 
Diffusiophoresis 
Electrostatic  effects 

-  Attraction 

-  Repulsion 

Gravitational  settling 

Hygroscopicity 

Impaction 

Momentum 

Physical  properties 

Resusoension 

Shape 

Size 

Solubility 
Thermophores is 

Gases 

Chemical  activity 
Diffusion: 

-  Brownian 

-  Eddy 

Partial  pressure  in 
equilibrium  with 
surface 
.  Solubility 


Accommodation 

-  Exudares 

-  Tfichomes 

-  Pubescence 

-  Wax 

Biotic  surfaces 
Canopy  growth: 

-  Dormant 

-  Expanding 
Senescent 
Canopy  structure: 

•  Areal  density 

-  Bark 

-  Bole 

-  Leaves 

-  Porosity 

-  Reproductive 

structure 

-  Soils 

-  Stem 

-  Type 

Electrostatic  properties 
Leaf -Vegetation; 

-  Boundary  layer 

-  Change  at  high  winds 

-  Flutter 

-  Stomatal  resistance 
Non-biotic  surfaces 
pH  effects  on; 

-  Reaction 

-  Solubility 
Pollutant  penetration 

and  distribution  in 
canopy 

Prior  deposition  loading 
Water 


*from  Sehmel,  1980 
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Figure  3.  Observed  deposition  velocities  as  a  function  of  particle  size  for 
1.5  g/cm^  density  particles.  Measured  by  Sehmel  and  Sutter  (1974) 
and  Holler  and  Schumann  (1970).  Figure  from  Sllnn  et  al.  (1978). 

Although  it  is  not  practical,  or  even  understood  how,  to  include  all  of  the 
variables  listed  in  table  1  into  the  deposition  model,  it  is  possible  to 
parameterize  many  of  the  most  important  effects  known  to  control  deposition 
rates,  using  only  directly-available  or  easily-computed  variables  which  describe 
the  state  of  the  atmosphere,  surface  conditions,  and  pollutant  properties. 

Many  models  of  dry  deposition  express  the  deposition  velocity  as  the  inverse  of 
a  sum  of  "resistances'*  plus,  for  particles,  gravitational  settling  velocity 
terms.  The  resistances  represent  the  opposition  to  transport  of  the  depositing 
material  from  a  reference  height  through  the  turbulent  atmospheric  surface  layer, 
and  through  a  quasi-laminar  layer  just  above  the  surface  to  the  surface  itself. 
The  major  processes  that  determine  these  resistances  are  briefly  described  below. 

3.1.1  Gravitational  settling 

The  gravitational  settling  velocity  is  a  function  of  the  particle  size,  shape, 
and  density.  Figure  4  shows  the  gravitational  settling  velocity  (V.j.  in  tfie 
figure)  as  a  function  of  particle  size  for  several  values  of  the  particle 
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density.  Note  that  the  gravitational  settling  velocity  represents  a  lower  limit 
to  the  deposition  velocity.  It  can  be  seen  that  for  larger  particles,  in  the 
range  of  20  to  40  in  diameter  and  higher,  the  deposition  velocity  approaches 
the  settling  velocity,  which  indicates  that  the  rate  of  deposition  is  dominated 
by  the  gravitational  settling  mechanism.  The  gravitational  settling  velocity 
decreases  with  decreasing  particle  size  and  density.  However,  for  particles 
smaller  than  about  20  pm  in  diameter,  the  deposition  velocity  curve  shows  larger 
and  larger  deviations  from  the  gravitational  settling  curve  as  the  particle  size 
decreases.  This  is  due  to  the  effect  of  other  mechanisms,  discussed  below,  in 
enhancing  the  deposition  rates  of  smaller  particles  above  those  predicted  by 
gravitational  settling  alone. 


Figure  4.  Predicted  deposition  velocities  for  u«  -  100  cm/s  and  particle 
densities  of  1,  4,  and  11.5  g/cm^.  Also  shown  is  the  gravitational 
settling  velocity  (V^) .  Fignre  from  Sehmel  (1980) . 
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3.1.2  Atmospheric  diffusion 


The  rate  of  deposition  can  sometimes  be  limited  by  the  transfer  of  pollutant 
material  to  the  vicinity  of  the  surface  by  atmospheric  turbulence.  For  example, 
this  would  typically  occur  during  very  stable  conditions  for  an  elevated  plume 
of  material  composed  of  small-sized  particles  with  small  gravitational  settling 
velocities.  In  the  lowest  layer  of  the  atmosphere,  the  aerodynamic  resistance 
is  used  to  parameterize  the  rate  of  mixing  in  terms  of  the  wlndspeed,  atmospheric 
stability,  and  surface  roughness  length.  The  aerodynamic  resistance  generally 
decreases  (that  is,  the  deposition  velocity  will  Increase)  with  increasing 
windspeed  and/or  surface  roughness. 

3.1.3  Quasi- laminar  layer. 

Over  smooth  surfaces,  a  thin  nonturbulent  sublayer  develops  that  can  be  a 
significant  obstacle  to  the  transfer  of  the  pollutant  onto  the  surface.  For 
rough,  real-world  surfaces,  this  sublayer  is  constantly  changing  and  is  likely 
to  be  intermittently  turbulent.  For  this  reason,  Hicks  (1982)  calls  this  layer 
the  "quasi -laminar"  layer.  It  is  also  know  as  the  deposition  layer.  Small 
particles  (<  0.05  fi,m  in  diameter)  are  transported  through  the  quas i -  laminar  layer 
primarily  by  Brownian  diffusion.  This  process  becomes  less  efficient  as  the 
particle  size  increases.  Particles  in  the  2-  to  20-A4m  diameter  range  tend  to 
penetrate  the  quasi-laminar  layer  by  inertial  impaction.  However,  particles 
larger  than  20  pm  in  diameter  are  less  efficiently  captured,  so  the  inertial 
impaction  mechanism  is  most  effective  in  the  2-  to  20 -pm  diameter  size  range. 
Because  particles  in  the  0.1-  to  1.0-pm  diameter  size  range  are  not  efficiently 
transported  across  the  quasi-laminar  layer  by  either  Brownian  diffusion  or 
inertial  impaction,  particles  in  this  size  range  have  the  lowest  deposition 
velocities . 

Several  candidate  deposition  models  that  parameterize  some  or  all  of  the  physical 
processes  described  above  have  been  selected  for  further  study.  Individual 
components  of  other  models  have  also  been  reviewed.  Ultimately,  it  may  be 
advantageous  to  combine  the  best  features  of  several  algorithms  into  a  hybrid 
model.  Brief  descriptions  of  the  various  models  for  computing  deposition 
velocities  of  particulate  matter  are  provided  below.  The  data  requirements  of 
all  the  models  considered  are  relatively  simple  (particle  size,  density,  surface 
roughness,  and  routine  meteorological  parameters  to  compute  the  friction  velocity 
and  Monin-Obukhov  length) .  Some  initial  results  of  sensitivity  testing  and 
intercomparisons  of  the  models  are  also  sununarized  below. 

3. 1.3.1  Sehmel  model.  Sehmel  (1980)  and  Sehmel  and  Hodgson  (1978)  proposed  a 
model  for  predicting  deposition  velocities  of  particles  above  smooth  surfaces. 
The  basis  of  the  model  is  a  set  of  wind  tunnel  observations  of  deposition  for 
monodispersed  particles  to  surfaces  such  as  gravel,  artificial  grass,  brass  shim 
stock,  and  water.  The  model  consists  of  empirical  equations  for  transfer 
resistances  derived  from  a  least-squared  empirical  fit  of  deposition  velocity  as 
a  function  of  particle  size,  density,  surface  roughness,  and  friction  velocity. 
The  equations  were  converted  into  a  computer  code  by  B.  Croes  of  the  California 
Air  Pesources  Board  (CARS)  and  is  sometimes  known  as  the  CARB  model.  It  is  also 
used  in  the  fugitive  dust  model  (FDM) . 
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In  the  Sehmel  approach,  integrated  resistances  to  mass  transfer  are  computed 
within  two  layers.  The  first  layer  extends  from  a  reference  height  of  1  m  down 
to  1  cm  above  the  surface.  In  this  layer,  atmospheric  turbulence  dominates  mass 
transfer.  Eddy  diffusivities  are  used  to  describe  the  transfer  rate.  The  second 
layer  is  the  deposition  surface  layer  within  one  centimeter  of  the  surface.  The 
integrated  resistance  within  the  deposition  surface  layer  is  derived  from  a 
statistical  fit  of  the  wind  tunnel  particle  deposition  data.  Sehmel  and  Hodgson 
express  the  deposition  velocity  as; 


V.  = 


-  exp[-Vy(X 


12 


(52) 


where 

V(j  is  the  deposition  velocity  (cm/s)  , 

Vg  is  the  gravitational  settling  velocity  (cm/s), 

1^2  is  the  atmospheric  diffusional  resistance  (dimensionless) , 

I3  is  the  surface  resistance  integral  (dimensionless),  and, 
u*  is  the  surface  friction  velocity. 

The  atmospheric  diffusion  resistance  used  by  Sehmel  and  Hodgson  (1978)  is  based 
on  the  flux  profile  relationships  of  Businger  et  al.  (1971).  For  neutral  or 
stable  conditions, 

Jj^2  *  z^)  +  -  Z2)/L)/k  (53) 

where 

Zi  is  the  upper  limit  of  the  atmospheric  diffusional  resistance  integral 
(that  is,  100  cm), 

Z2  is  the  lower  limit  of  integral  (that  is,  1  cm), 

L  is  the  Monin-Obukhov  length  (cm),  and, 

k  is  the  von  Karman  constant  (Sehmel  used  a  value  of  0.35). 

For  unstable  conditions,  the  atmospheric  diffusional  resistance  integration 
yields : 


(Inl  +  i-)J  ) 

2  .  [tan'^(A^2^  “  tan’^Cfl^pi/Jc 

(54) 

=  (1,  -  15.Zi/l)°-25 

(55) 

=  (i-  -  is.z^/l'P-^^ 

(56) 
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The  surface  resistance  integral  is  an  empirical  relationship  based  on  wind  tunnel 
observations.  For  particles  with  a  diameter  >0.01  mm, 

J  =  exp  {-378.051  +  I6.4981n(5c)  + 
ln{t)  [  -11.818  -  0.2863  ln(t  )  + 

0.3226ln(d  10  /z  )  -  0.3385  ln(I?/z  u  )  )  ]  - 
12,804ln(dl0  )} 


where , 

D  is  the  Brownian  diffusion  coefficient  (cm  /s) , 

t  is  a  relaxation  time  (dimensionless) , 

Sc  is  the  Schmidt  number, 

z  is  the  surface  roughness  length  (cm) ,  and, 

d  is  the  particle  diameter  (mm) . 

Sehmel's  formulation  generally  shows  a  reasonable  variation  of  deposition 
velocity  as  a  function  of  its  model  parameters  (that  is,  density,  size,  surface 
roughness,  friction  velocity) .  The  predicted  deposition  velocity  is  close  to  the 
gravitational  settling  velocity  for  large  particles  (for  example,  greater  than 
about  20  mm  diameter) ,  and  decreases  with  decreasing  particle  size  to  about 
0,1  mm  to  1.0  mm,  where  it  reached  a  minimum.  The  deposition  velocity  then 
increases  with  decreasing  particle  size  for  smaller  sized  particles.  This 
behavior  is  consistent  with  the  importance  of  Brownian  motion  in  enhancing 
deposition  rates  for  very  small  particles.  The  Sehmel  scheme  produces  increased 
deposition  rates  for  increased  particle  density,  surface  roughness  lengths,  and 
friction  velocities.  Recent  laboratory  studies  by  Noll  and  Fang  (1989)  suggest, 
however,  that  Sehmel's  model  may  underpredict  deposition  velocities  fo  particle 
diameters  greater  than  about  30  mm. 

A  main  limitation  of  the  scheme  stems  from  the  lack  of  generality  of  the  highly 
empirical  relationship  for  the  surface  resistance  integral  (equation  (57)).  I 
is  based  on  wind  tunnel  data  for  relatively  smooth  surfaces  under  a  limited  range 
of  conditions.  For  example,  in  order  to  avoid  extrapolation,  the  GARB  implemen¬ 
tation  of  the  model  does  not  allow  the  surface  roughness  length  used  in  the 
algorithm  to  exceed  10  cm,  even  though  many  real  world  surfaces  have  signifi¬ 
cantly  greater  roughness  lengths.  In  addition,  sensitivity  testing  of  the  model 
has  shown  that  it  exhibits  some  nonphysical  behavior  when  the  inputs  are  varied 
beyond  the  range  of  conditions  tested  in  the  wind  tunnel.  For  example,  the 
Sehmel  model  shows  a  very  strong  sensitivity  to  temperature  which  is  not 
exhibited  by  other  deposition  models.  For  particles  in  the  0.1-  to  1.0-Mm 
diameter  size,  a  change  of  nearly  an  order  of  magnitude  in  deposition  velocity 
was  predicted  for  a  temperature  change  from  0  “F  to  100  “F  (figure  5).  This 
nonphysical  behavior  is  probably  an  artifact  of  the  regression  equations  used  to 
fit  the  surface  resistance  integral  to  the  wind  tunnel  data.  The  original  wind 
tunnel  tests  were  performed  at  a  constant  temperature,  and,  as  a  result,  applica¬ 
tion  to  a  realistic  range  of  atmospheric  conditions  involves  extrapolation 
outside  the  range  on  which  the  model  was  developed. 
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A  second  problem  noted  from  the  sensitivity  testing  is  chat  the  deposition 
velocity  curves  display  a  kink  at  a  particle  diameter  of  about  0.03  pim.  For 
particles  smaller  than  about  0.1  fiW  in  diameter,  the  deposition  velocity 
increases  with  decreasing  particle  diameter.  However,  the  Sehmel  model  shows 
this  trend  only  to  about  0.03  pm  in  diameter,  beyond  which  the  deposition 
velocity  is  predicted  to  decrease  or  level  off  with  decreasing  particle  diameter 
(see  figure  6) .  These  kinks  in  the  deposition  velocity  curves  appear  to  be 
another  artifact  of  use  of  the  regression  equations  outside  the  range  of 
conditions  on  which  the  model  was  developed.  These  limitations  have  some 
implications  on  the  ability  of  this  model  to  handle  the  entire  range  of  likely 
applications . 


PARTICLE  DIAMETER  (MICRONS) 

Figure  6.  Deposition  velocity  as  a  function  of  particle  diameter  as  predicted 
by  the  Sehmel/CARB  model  for  surface  roughness  lengths  of  0.001, 
0.01,  3,  and  10  cm.  The  particles  have  an  assumed  density  of  1  g/cc 
and  are  in  a  neutral  flow  with  u*  -0.1  m/s . 
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3. 1.3. 2  Resistance  models.  Resistance -based  models  for  particle  deposition  have 
been  used  in  many  models,  including  ADOM  (Pleim  et  al.,  1984),  GALPUFF  (Scire  et 
al.,  1990),  CALGRID  (Yamartino  et  al.,  1989),  and  a  new  version  of  the  Urban 
Airshed  Model,  UAM-V.  The  basic  approach  used  in  these  models  is  similar,  with 
special  enhancements  or  features  designed  to  accommodate  the  specific  focus  of 
the  model  (for  example,  oxidants,  secondary  particulate  matter,  toxic 
pollutants) .  Although  some  of  these  models  are  mesoscale  or  regional  scale 
Eulerian  grid  models  (for  example,  ADOM,  CALGRID),  their  deposition  formulations 
are  suitable  for  application  to  Gaussian  plume  or  Lagrangian  particle  models  on 
much  smaller  horizontal  spatial  scales. 


The  general  approach  used  in  the  resistance  models  is  to  include  explicit 
parameterizations  of  the  effects  of  Brownian  motion,  inertial  impaction,  and 
gravitational  settling.  The  deposition  velocity  is  written  as  the  inverse  of  a 
sum  of  resistances  to  pollutant  transfer  through  the  atmosphere  to  the  surface, 
plus  gravitational  settling  terms  (Slinn  and  Slinn,  1980;  Pleim  et  al.,  1984). 


= 


(58) 


where 

Vj  is  the  deposition  velocity  (cm/s), 

Vg  is  the  gravitational  settling  velocity  (cm/s) , 

V,  is  the  aerodynamic  resistance  (s/cm) ,  and 

Vd  is  the  deposition  layer  resistance  (s/cm). 

Note  that  for  large  settling  velocities,  the  deposition  velocity  approaches  the 
settling  velocity  (v^  -*  L  Vg)  ,  whereas,  for  small  settling  velocities,  Vj  tends 
to  be  dominated  by  the  r,  and  r^  resistance  terms. 

The  lowest  few  meters  of  the  atmosphere  can  be  divided  into  two  layers  --  i  fully 
turbulent  region  where  vertical  fluxes  are  nearly  constant,  and  the  thin  quasi - 
laminar  sublayer.  The  resistance  to  transport  through  the  turbulent,  constant 
flux  layer  is  referred  to  as  the  aerodynamic  resistance.  It  is  usually  assumed 
that  the  eddy  diffusivity  for  mass  transfer  within  this  layer  is  similar  to  that 
for  heat.  The  atmospheric  resistance  formulation  used  in  all  of  the  resistance 
models  is  based  on  Wesely  and  Hicks  (1977). 


where 

is  a  stability  adjustment  factor, 
u*  is  the  surface  friction  velocity, 
k  is  the  von  Karman  constant  (0.4), 
z  is  ’-.eight  (m) ,  and 
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Zg  is  the  surface  roughness  length  (m) . 

The  parameterization  of  the  deposition  layer  resistance  terms  varies  somevrhat. 
The  approach  used  by  Pleim  et  al.  (1984)  is: 

=  [  5C-2/3  +  <60) 


where 

Sc  is  the  Schmidt  number  (Sc  -  I'/D) , 

1/  is  the  viscosity  of  air, 

D  is  the  Brownian  diffusivity  of  the  pollutant  in  air, 

St  is  the  Stokes  number  [St  -  {v^/ g)  (u^/v)  ],  and 

g  is  the  acceleration  due  to  gravity  (9.8  m/s^) . 

The  first  term  of  equation  (60),  involving  the  Schmidt  number,  parameterizes  the 
effects  of  Brownian  motion.  This  term  controls  the  deposition  rate  for  small 
particles.  The  second  term.  Involving  the  Stokes  number,  is  a  measure  of  the 
imporiance  of  inertial  impaction.  The  inertial  impaction  term  tends  to  dominate 
for  intermediate-sized  particles  in  the  2-  to  20-/im  diameter  size  range. 

The  deposition  formulation  in  UAM-V  is  still  under  development,  and  is  currently 
being  reviewed  by  Sigma  Research  under  the  Lake  Michigan  Oxidant  Study  (LMOS) 
program.  Its  main  difference  from  that  of  Pleim  et  al.  (1984)  lies  in  the 
formulation  of  the  inertial  impaction  term  of  r^.  That  is, 

Tj  =  [Ci  Sc‘2''^  +  Cj  5t]‘^  uV' 

where  Cj,  C2  are  constants,  linearizes  the  appearance  of  the  Stokes  niimber. 

These  models  produce  a  pattern  of  deposition  velocity  as  a  function  of  particle 
size  and  density,  similar  to  that  of  Sehmel.  However,  the  ADOM/CALPUFF/  CALGRID 
approach  of  Pleim  et  al.  ,  tends  to  predict  somewhat  higher  deposition  velocities 
in  the  5-  to  15 -pm  diameter  size  range  than  the  Sehmel  model  and  lower  values  in 
the  0.1  to  5.0-pm  diameter  range.  Although  the  general  shape  of  the  deposition 
velocity  curves  are  similar,  the  minimum  deposition  velocity  in  the  Pleim  et  al . 
formulation  tends  to  occur  at  larger  particle  diameters  (that  is,  closer  to  L.O 
pm  than  0,1  pm)  than  in  the  Sehmel  model.  The  constants  and  C2  in  the  UAM-V 
scheme  were  derived  to  force  a  minimum  deposition  velocity  at  0.2  pm  diameter. 

The  parameterizations  of  the  resistance  models  for  Brownian  motion  and  inertial 
impaction  effects  involve  empirical  factors  derived  from  a  number  of  field  and 
wind  tunnel  studies.  The  sensitivity  analyses  with  these  models  shows  no  unusual 
rasponse  to  temperature  variations  over  the  ambient  temperature  range  from  0  T 
to  100  °F.  These  resistance  models  also  all  show  a  steady  increase  of  deposition 
velocity  with  decreasing  particle  diameter  in  the  small  particle  range  (<  0,1  pm 
in  diameter) . 
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Figures  7  thru  10  display  predictions  of  the  UAM  deposition  velocity  model  for 
particle  sizes  ranging  from  10"^  fim  to  20  fim  and  for  seven  different  values  of 
surface  roughness,  z^,  ranging  from  10"^  cm  to  100  cm.  Figures  7  and  8  are  for 
stable  meteorological  conditions  with  friction  velocity  values  of  0.1  m/s  and  0.5 
m/s,  respectively,  whereas  figures  9  and  10  involve  the  assumption  of  convective 
conditions.  As  most  of  the  deposition  variability  is  due  to  u*,  the  influences 
of  surface  roughness  and  stability  per  se  are  seen  to  be  relatively  small. 


PARTICLE  DIAMETER  (MICRONS) 


Figure  7.  Deposition  velocity  versus  particle  size  as  predicted  by  the  UAM-V 
model  for  surface  roughnesses  of  10"’,  10"^,  3,  10,  20,  50,  and  100 
cm.  The  model  predicted  deposition  velocities  are  for  particles  with 
a  density  of  1  g/cc  in  a  flow  with  stable  stratification  and  a 
friction  velocity  of  0.1  m/s. 
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PEPOSITION  VELOCITY  ICM/SECI 


PARTICLE  DIAMETER  (MICRONS) 


Figure  8.  Deposition  velocity  versus  particle  size  as  predicted  by  the  UAM-V 
model  for  surface  roughnesses  of  10“^,  10“^,  3,  10,  20,  50,  and  100 
cm.  The  model  predicted  deposition  velocities  are  for  particles  with 
a  density  of  1  g/cc  in  a  flow  with  stable  stratification  and  a 
friction  velocity  of  0.5  m/s. 
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DEPOSITION  VELOCITY  I  CM/SEC) 


Figure  9 .  Deposition  velocity  versus  particle  size  as  predicted  by  the  UAM-V 
model  for  surface  roughnesses  of  10“^,  10"^,  3,  10,  20,  50,  and  100 
cm.  The  model  predicted  deposition  velocities  are  for  particles  with 
a  density  of  1  g/cc  in  a  flow  with  unstable  stratification  and  a 
friction  velocity  of  0.1  m/s. 
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DEPOSITION  VELOCITY  (CH/SEC) 


Figure  10.  Deposition  velocity  versus  particle  size  as  predicted  by  the  UAM-V 
model  for  surface  roughnesses  of  10"^,  10“*^,  3,  10,  20,  50,  and  100 
cm.  The  model  predicted  deposition  velocities  are  for  particles  with 
a  density  of  1  g/cc  In  a  flow  with  unstable  stratification  and  a 
friction  velocity  of  0.5  m/s. 


3.2  Deposition  Removal  in  Particle  Models 

There  are  two  widely  used  methodologies  for  including  dry  deposition  of 
particles  - -discrete  particle  removal  and  probabilistic  reduction  of  particle 
mass.  Both  methods  involve  the  computation  of  capture  probabilities  for  the 
surface  to  ensure  that  the  deposition  flux,  F,  matches  that  predicted  by  the 
relation,  F  -  v^  •  CCz^gf)  where  C(Zr,f)  is  the  concentration  at  the  same  reference 
height  used  in  the  computation  of  the  deposition  velocity,  v<j. 

In  the  case  of  discrete  particle  removal,  particles  with  a  trajectory  time  step 
end  point  below  the  surface,  that  is,  z,  <  0  are  either  reflected,  by  reassigning 
this  coordinate  as  -z,,  or  absorbed  to  provide  an  appropriate  deposition  flux. 
Overcamp  (1976)  introduced  a  partial  reflection  formulation  for  Gaussian  plume 
models  including  both  dry  deposition  and  gravitational  settling,  while  Rao  (1981) 
has  developed  a  formulation  based  on  an  analytic  solution  of  the  equivalent 
gradient  transfer  theory  problem.  Janicke  (1985)  has  extended  the  notion  of 
partial  reflection  to  Lagrangian  particle  models  but  has  recently  concluded 
(Janicke,  1990)  that  the  notion  of  applying  weights  of  zero  or  one  to  particle 
survival  probabilities  serves  to  increase  the  statistical  uncertainty  associated 
with  counting  particles  (that  is,  counting  statistics  errors)  for  evaluating 
concentrations.  In  his  LASAT-C  model  (Janicke,  1990),  all  particles  which  reach 
the  z  -  0  surface  are  reflected  and  a  fractional  survival  probability,  p,,  (0  < 
Ps  <  1)  is  applied  to  all  such  particles.  The  probabilities  so  determined 
involve  several  velocity,  time,  and  distance  scales  and  the  formalism  is  not 
obviously  transferable  to  the  envisioned  KS  turbulence  fomnulation.  In  addition, 
such  a  scheme  may  continue  to  distribute  capture  probabilities  across  relatively 
few  particles  with  higher  attendant  computational  noise. 

Zannetti  and  Al-Madani  (1983)  proposed  a  probabilistic  method  which  also  is 
applied  only  to  those  particles  which  have  passed  below  the  z  -  0  surface.  In 
their  model,  the  fraction  of  particles  which  are  deposited  is  specified  as  p^, 
where 


-  1  -  exp{-At/T^)  (62) 

and  Tj  is  the  deposition  time  scale  which  is  a  function  of  all  the  variables 
influencing  deposition.  An  analogous  method  is  also  developed  to  describe 
particle  resuspension.  Such  a  method  is  interesting  both  because  it  treats 
deposition  as  an  exponential  process  and  because  it  can  be  implemented  either  in 
discrete  particle  removal  schemes  or  in  mass  reduction  schemes.  In  probabilisti 
mass  reduction  schemes,  the  particle  is  not  removed  but  instead  assigned  a 
reduced  mass , 

M{C  +  At)  =  (1  -  Pj)-Af(t)  ,  (63) 

for  use  in  subsequent  concentration  calculations.  This  representation  is  more 
explicit  in  indicating  that  the  fractional  mass,  p*,  (or  probability  of) 
surviving  the  deposition  process  is  simply: 

Ps  =  1  -  Pd  =  expC-Ac/Tj)  ,  (6^) 

and  therefore  treatable  as  an  exponential  process. 
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Equation  (64)  is  also  the  form  used  to  describe  the  surviving  fraction  of 
material  in  one -layer,  long  range  transport  models  or  in  the  lowest  layer  of 
multi-layer  grid  models  and  results  directly  from  the  solution  of  the  first-order 
differential  equation  for  material  loss: 


dC 

dt 


=  -F-A 


(65) 


where 


V  is  the  volume  of  the  cell, 

A  is  the  cell's  area  in  the  x  -  y  plane, 

and 

F  is  the  deposition  flux  given  as 

F  -  Vj«C  in  terms  of  the  deposition  velocity  Vj. 

This  differential  equation  is  easily  solved  to  yield, 

p,  =  C{At) /C(0)  »  exp(-Vj-Ac///)  ,  (66) 

where  At  is  the  duration  of  the  time  step,  and  H  -  V/A  is  the  cell's  vertical 
depth  over  which  the  pollutant  is  assumed  well  mixed. 

Equation  (66)  raises  the  interesting  question  of  why  should  an  arbitrarily 
defined  cell  depth  H  affect  the  deposition  of  material?  In  fact,  it  is  easy  to 
show  that  provided  the  cell  is  not  excessively  depleted  (that  is,  Vj,At/H  «  1) 
the  number  of  particles  deposited,  n^,  from  a  cell  containing  n^  total  particles 
is  just 

~  Hj.  '  iv^'Lt/H)  ;  (67) 

however,  the  number  of  particles  in  the  cell  equals  the  number  density 
concentration.  C/m  ,  (where  mg  is  mass  per  particle)  times  the  cell  volume,  AWH, 
to  yield 

-  (C/mJ  'A'v^'^t  =  {F/nio)  'A'^t  (68) 

thus,  equation  (68)  shows  that  for  <}»  ■  v^-AT/H  small  enough  that  1-exp  -  4* 
the  depth  of  the  cell  is  not  of  importance. 

We  now  extend  this  concept  to  develop  a  model  for  particle  survival  probabi  lity 
that  is  a  reasonable  function  of  height,  z,  above  the  surface.  The  simplest  such 
model  that  preserves  the  cumulative  exponential  nature  of  the  process  is 

Pgiz)  =  A-exp(=P  Vjj-At/z)  ,  (69) 
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where  A  is  an  overall  normalization  of  order  unity  and  ^  is  a  constant 
coefficient.  The  primary  constraint  that  equation  (69)  must  satisfy  is  that  its 
average  value  over  the  interval  0  s  z  <  H  corresponds  with  the  value  of  p,  from 
equation  (66) . 

Evaluation  of  this  averaging  integral  yields 

Pa  =  %  ^  dz  Pg  iz)  =A£'2(P<|))  <705 

n  o 

where  Ez  is  the  exponential  integral  function  defined  as 

{X)  =  Idt  <71) 

1 


and  numerical  values  of  0  in  the  range  of  to  i  leads  to  the  desired 
approximate  equality 

E2  (P<t>)  -  e‘*  for  4>  <  1.  <72) 

Thus,  to  make  the  equation  (69)  formulation  most  generally  applicable,  we  define 
the  overall  normalization  as 

A  3  e-<|)  £’2<P‘I>) 


where,  as  before,  <|>  *  Examination  of  equation  (73)  using  expansions  of 
the  E2  function  and  for  appropriate  values  of  0  (for  example,  P  -  1/6)  show  that 
A  is  always  less  than  unity  so  that  one  is  never  in  the  embarrassing  and  puzzling 
position  of  having  p8(z)  values  greater  than  one. 


Another  advantage  of  the  equation  (69)  formulation  is  that  the  cumulative 
survival  probability,  p^,  after  n  time  steps  can  be  written  as 


n 


Pg  =  n  A  exp  (-pVj,-Ar/Z_i) 
i=l 


=A  exp 


n 


E  (l/Zi) 

i=l 


(74) 


Thus,  for  computation  of  a  particle's  cumulative  survival  probability  it  becomes 
only  necessary  to  accumulate  the  harmonic  sum  of  individual  particle 
trajectory  values,  or  possibly  vj^/z^  values  if  the  deposition  velocity  varies 
sufficiently  rapidly  in  space  or  time.  Actual  evaluation  of  equation  (74) 
further  demands  that  values  of  z^  >  H  be  excluded  from  the  stim  and  that  Zj  -  0 
values  be  replaced  with  small  positive  values,  perhaps  of  order  10"^  •  H  for  a 
typical  H  “  20  m . 

A  further  advantage  of  the  proposed  formulation  is  that  there  are  no  remaining 
unspecified  time,  length,  or  velocity  scales  to  contend  with.  Finally,  we  note 
that  since  mathematical  particles  are  never  eliminated,  the  same  entity  can 
simultaneously  represent  physical  particles  from  a  range  of  particle  size 
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categories  (not  strongly  influenced  by  gravitational  settling) .  It  is  only 
necessary  to  track  values  of  survival  probability  separately  for  different  v^ 
values . 


3.3  Point  and  Path- Integrated  Concentration  Computations 


Computation  of  instantaneous  point  and  path-averaged  concentrations  requires  the 
"counting"  of  particles  within  a  volume  region  about  the  point  or  path  of 
interest.  These  volumetric  sampling  regions  may  either  be  defined  in  a  "box  car" 
sense,  with  particles  within  the  box  being  assigned  a  weight  of  one  and  chose 
outside  assigned  a  weight  of  zero,  or  in  a  diffuse  volume  sense  where  particles 
are  assigned  a  continutjtm  of  weights  depending  on  their  proximity  to  the  actual 
sampling  location.  A  simple  example  of  Che  latter  concept  involves  the  second- 
dimensional  or  third- dimensional  Gaussian  weighting  function 

w(z)  a  — -  *  exp  ,  (75) 

a)  12  / 

where  F  is  the  scalar  distance  between  the  particle  and  sampling  line  (n-2)  or 
point  (n-3)  and  >J2no  is  the  equivalent  box  normalization  length  scale.  In  the 
case  of  the  line  integral,  the  length  of  the  finite  line,  L,  will  serve  as  the 
third  dimensions  for  defining  the  volxune  of  Che  cylindrical  region.  For 
L  i  V2iio  the  volume  associated  with  the  line  integration  will  exceed  chat 
associated  with  the  point  sample,  leading  typically  to  a  greater  number  of 
sampled  particles  and  consequently  a  lower  statistical  uncertainty  associated 
purely  with  sample  size,  N.  This  statistical  uncertainty,  o^,  is  simply  obtained 
from  the  theory  of  Poisson  statistics  and  yields  the  fractional  uncertainty 

l/y/N  for  >  1 
N  ^ 


for  particles  with  uniform  weights  w^.  For  a  mixture  of  weights,  the  effective 
sample  size  is  smaller  than  N  and  is  given  as 


N 


|l/2 


fox  N  ■>  1 


(77) 


The  implication  of  equation  (7b)  is  that  to  reduce  the  statistical  counting 
uncertainty  to  say  10  percent,  one  must  count  100  particles  (or  more  if  w^  vary) 
within  the  sampling  volume.  These  considerations  place  a  strong  constraint  on 
the  total  number  of  particles  which  must  be  released  to  achieve  a  given  level  of 
uncertainty.  Attempts  to  reduce  computation  time  by  limiting  the  number  of 
particles  lead  one  to  increase  the  sampling  volume  length  scales,  s  and  L; 
however,  increasing  these  scales  is  tantamount  to  reducing  the  spatial  resolution 
of  the  model. 

The  above  methodology  is  very  nearly  equivalent  to  the  kernel  sampling  methods 
developed  by  Lorimer  (1986)  and  others,  except  that  the  spatially  diffuse 
weighting  function  (for  example,  that  of  equation  (75)  for  n-3)  is  now  attached 
to  each  moving  particle  rather  than  the  sampler.  In  kernel  theory  the  length 
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scale,  a,  is  allowed  to  grow  in  time  with  a  growth  rate  and  law  which  may 
encompass  varying  fractions  of  the  turbulent  eddy  spectrxun.  This  is  tantamount 
to  ensemble  averaging  over  possible  realizations  for  that  fraction  of  the 
turbulent  spectra  included  and  may  create  difficulty  in  interpreting  results 
especially  when  a  hybrid  turbulent  parameterization  (for  example,  KS  theory  plus 
the  Langevin  equation)  is  implemented.  Yamada  and  Bunker's  (1988)  kernel 
estimator  in  the  RAPTAD  model  utilizes  the  full  spatial  a  values  from  Taylor's 
(1921)  homogeneous  theory  and  therefore  will  yield  results  as  highly  smoothed  as 
those  resulting  from  a  Lagranglan  puff  model. 

4.  DEVELOPMENT  AND  TESTING  OF  A  PROTOTYPE  MODULE 

4.1  Driver  Program  and  Mean  Flow  Trajectory  Module 

Referring  back  to  figure  1  it  is  evident  chat  a  main  driver  program  is  needed  to; 

•  read  model  control  files ; 

•  access  meteorological  data  files  and  additional  external  data  files; 

•  serve  as  an  envelope  for  modules  for 

source  characterization  and  particle  release, 

synthetic  KS  turbulence  field  creation,  temporal  evolution  and 
sampling, 

-  particle  trajectory  simulation, 
particle  deposition  and  removal, 

point  and  path  integrated  instantaneous  and  averaged 
concentrations ,  and 

-  multispectral  transmission  attenuation;  and 

•  produce  appropriate  statistical  summary  and  graphical  output  files. 

A  basic  driver  program  has  been  developed  and  tested.  The  present  version  reads 
the  basic  model  control  files,  reads  gridded  hourly  fields  of  mean  winds  (that 
is,  third- dimensional  fields  of  u,  v,  w)  and  other  meteorological  variables  (for 
example,  second- dimensional  fields  of  mixing  depth,  u*,  w*,  Monin-Obukhov  length) 
generated  by  the  CALMET  meteorological  model,  controls  the  release  of  particles 
from  a  user-specified  source,  and  generates,  stores,  updates,  and  outputs 
trajectories  of  particles  driven  by  the  mean  flow.  This  driver  program  also 
provides  a  convenient  platform  for  insertion  of  the  above  mentioned  modules- -some 
of  which  have  been  developed  and  tested  separately. 

Figure  11  displays  the  trajectories  of  five  particles  released  600  s  apart  into 
a  CALMET  generated  surface  layer.  The  particles  are  tracked  for  5  h  with  arrows 
indicating  particle  position  and  direction  every  300  s.  In  this  example, 
particles  respond  solely  to  a  local  interpolated  mean  flow  determined  via 
bilinear  interpolation  in  x,  y,  and  linear  interpolation  in  z  between  the 
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adjacent  vertical  levels.  The  interpolation  procedures  and  other  discrete 
operators  have,  of  course,  been  placed  into  individual  FORTRAN  77  subroutines  and 
functions  to  yield  a  readable,  structure  code. 


Y 


X 


Figure  11.  Trajectories  of  five  particles  released  600  s  apart  into  a  CALMET 
mesoscale  wind  field.  Arrows  depict  particle  position  and  direction 
at  300  s  Intervals.  The  x-y  region  shown  is  5  km  by  20  km. 
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The  graphical  output  postprocessor,  which  generated  the  x-y  plot  in  figure  11, 
can  also  plot  the  trajectories  in  other  planes.  Figures  12  and  13  show  the  same 
particle  trajectories  in  the  x-z  and  y-z  planes,  respectively.  As  previously 
noted,  particles  follow  only  the  interpolated  hourly  mean  wind  fields.  Thus,  all 
five  particles  begin  on  the  same  trajectory  but  jump  to  different  trajectories 
as  the  next  time  step's  wind  fields  are  utilized.  This  behavior  would  be  altered 
by  implementing  a  time  interpolation  scheme  between  multiple  time  levels  of 
winds,  but  the  net  effect  would  still  involve  an  effective  "dispersion"  of 
particles  due  to  the  particle  release  time  differential. 


X 


Figure  12.  Trajectories  of  five  particles  released  600  s  apart  into  a  CALMET 
mesoscale  wind  field.  Arrows  depict  particle  position  and  direction 
at  300  s  intervals.  The  x-z  region  shown  is  5  loa  by  20  m. 
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Figure  13.  Trajectories  of  five  particles  released  600  s  apart  into  a  CALMET 
mesoscale  wind  field.  Arrows  depict  particle  position  and  direction 
at  300  s  intervals.  The  y-z  region  shown  is  20  km  by  20  m. 

4.2  Development  and  Initial  Testing  of  the  Kinematic  Simulation  Turbulence 
Module 

Individual  modules  have  been  developed  and  tested  for  the  initialization,  time 
advancement,  and  velocity  field  sampling  of  the  sub-grid-scale,  KS  flow  fields 
described  in  section  2.4.  Completely  separate  KS  fields  can  be  generated  for 
each  of  the  mesoscale  model's  (x,  y,  z)  grid  cells  (that  is,  cells  denoted  by 
indices  i,  j,  k)  ;  however,  some  meteorological  situations  such  as  the  convective 
boundary  layer  might  he  more  realistically  simulated  by  a  single  set  of  KS  fields 
for  each  x,  y  (or  i,  j)  column  below  the  mixed  layer  height. 
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The  KS  modules  have  been  conveniently  split  into  the  three  components  of 
initialization,  time  advancement,  and  field  sampling  so  that: 

•  a  field  containing  N  spectral  components,  with  N  arbitrary,  and  user 
specified  spectral  shape  can  be  initialized  for  an  arbitrary  second- 
dimensional  and  third-dimensional  mesoscale  grid; 

•  the  state  of  these  fields  may  be  stochastically  time  advanced  at  any 
time  interval,  independent  of  particle  release  rate  or  trajectory 
Increment  time  steps,  with  phase  evolution  factors  (that  is,  the  and 
T„  of  equations  (46)  and  (47)  also  established  in  the  initialization 
algorithm;  and 

•  the  time  analytic  values  of  the  components,  (u'  ,  v'  ,  w')  may  be  sampled 

in  space  and  time  as  frequently  as  dictated  by  the  number  of  particles 

or  other  parameters  of  the  problem. 

Thus,  we  have  constructed  a  flexible  set  of  KS  velocity  field  operator  functions 
that  can  be  applied  under  a  wide  range  of  conditions  and  application  circum¬ 
stances.  To  date  only  the  equation  (44)  homogeneous  turbulence  expressions  for 

an  vmbounded  domain  have  been  included  into  the  algorithms,  but  transition  to  an 
alternative  set  of  equations,  such  as  equations  (49)  through  (51)  would  only 
require  replacing  about  20  lines  of  code.  A  "switch"  might  be  ultimately 
supplied  to  choose  between  several  KS  field  formulations . 

In  the  demonstration  tests  which  follow,  the  objective  was  to  fill  a  single 
column  of  the  mesoscale  model  (that  is,  with  Ax  -  Ay  -  5  km  and  a  boundary  layer 
depth  of  Az  2  km)  with  a  KS  field  determined  by  only  five  spectral  components. 
These  spectral  components,  chosen  to  have  equal  variance  and  with  k  values 
beginning  at  2ir/Ax  (or  2jr/Az)  and  increasing  by  factors  of  three  for  each 
component,  were  envisioned  as  leading  to  an  unrealistic,  multi-spiked  spectrum 
of  an  unphysical  nature.  A  typical  time  series  of  sampled  u’  values  every  3  s 
(for  example,  figure  14)  suggests  several  dominant  frequencies.  However,  the 
fixed-point  Eulerian  power  spectra  produced  during  this  one  hour  trial  is 
remarkably  smooth  beyond  the  lowest  driving  frequency  and  displays  a  spectral 
energy  falloff  exponent  of  nS(n)  in  the  inertial  subrange  of  somewhat  less  than 
1.0  (for  example,  -  -0.7  to  -0.8)  for  u'  and  w'  versus  theoretical  values  of 
-2/3).  These  two  spectra  for  u'  and  w'  are  displayed  as  figures  15  and  16, 
respectively . 

Correlograms  for  these  u'  and  w'  time  series  are  displayed  in  figures  17  and  18, 
respectively,  as  a  function  of  time  lag.  The  u'  and  w'  correlograms  are 
remarkably  similar,  with  Lagrangian  time  scales  (same  as  Eulerian  time  scales  as 
mean  winds  were  set  to  zero)  of  about  T^  -  60  s. 

Unfortunately,  the  u'  data  are  found  to  yield  a  somewhat  smaller  Tl  than  that  of 
the  w'  data;  however,  this  is  a  direct  consequence  of  the  equation  (44) 
formulation  of  the  wind  components  in  terms  of  a  cross  product  of  vectors  a,  b, 
and  k.  Thus,  the  u’  component  involves  the  y  and  z  components  of  amplitudes  a, 
b,  whereas  the  w'  velocity  involves  the  x  and  y  components,  and  attempts  to 
assign  the  shorter  Lagrangian  time  scale  to  a^,  b^  amplitudes  prove  obviously 
counterproductive.  This  suggests  further  that  a  more  direct  formulation  of  the 
amplitudes  (that  is,  not  in  terms  of  cross  products),  as  provided  by  equations 
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Figure  14.  Time  series  of  KS  predicted  u'  velocity  fluctuations.  The  KS 
amplitude  and  phase  components  were  advanced  every  3  s  and  the 
resulting  winds  sampled  at  the  same  time  interval. 
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Figure  15.  Frequency  scaled  power  spectra,  nS(n),  for  the  u'  component.  The 
data  sample  consisted  of  a  1-h  record  of  KS  generated  winds  with  KS 
field  time  advancement  and  sampling  intervals  of  3  s . 
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Figure  16.  Frequency  scaled  power  spectra,  nS(n),  for  the  w'  component.  The 
data  sample  consisted  of  a  1-h  record  of  KS  generated  winds  with  KS 
field  time  advancement  and  sampling  Intervals  of  3  s . 
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Figure  17 .  Autocorrelation  function  for  the  u'  velocity  component  versus  lag 
time.  The  data  sample  consisted  of  a  1-h  record  of  KS  generated 
winds  with  KS  field  time  advancement  and  sampling  intervals  of  3  s . 
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IS-  Autocorrelatilon  function  for  the  w'  velocity  component  versus  lag 
time.  The  data  sample  consisted  of  a  1-h  record  of  KS  generated 
winds  with  KS  field  time  advancement  and  sampling  intervals  of  3  s . 

Unfortunately,  the  u'  data  are  found  to  yield  a  somewhat  smaller  than  that  of 
the  w'  data;  however,  this  is  a  direct  consequence  of  the  equation  (44) 
formulation  of  the  wind  components  in  terms  of  a  cross  product  of  vectors  a,  b, 
and  k.  Thus,  the  u'  component  involves  the  y  and  z  components  of  amplitudes  a, 
b,  whereas  the  w'  velocity  involves  the  x  and  y  components,  and  attempts  to 
assign  the  shorter  Lagrangian  time  scale  to  a^,  b^  amplitudes  prove  obviously 
counterproductive.  This  suggests  further  that  a  more  direct  formulation  of  the 
amplitudes  (that  is,  not  in  terms  of  cross  products),  as  provided  by  equations 
(49)  thru  (51),  will  enable  greater  control  of  the  spectral  and  time  series 
characteristics  of  each  of  the  sub-grid-scale  velocity  components. 
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Finally,  Che  above  KS  winds  were  used  ».o  produce  particle  trajectories  and 
analyze  the  dispersion  produced  by  these  KS  winds  alone  or  in  combination  with 
higher  frequency  turbulence  components  characterized  by  a  Langevin  equation. 
Figure  19  shows  the  1-h  trajectory  of  a  single  particle  driven  purely  by  the  KS 
winds.  As  a  time  step  of  At  -  3  s  was  used  for  both  the  KS  field  update  and 
particle  transport,  particle  position  and  direction  arrows  are  not  shown  (that 
is,  as  they  are  not  individually  distinguishable). 


Figure  19.  Trajectory  of  a  single  particle  subjected  to  the  KS  generated 
sub-grid-scale  flow  fields.  The  particle  trajectory  was  advanced  in 
time  steps  of  3  s  for  1-h.  The  x-y  region  shown  is  1  km  by  1  km. 
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When  10  particles  were  released  3  s  apart  into  the  KS  winds,  interparticle 
separation  grew  rather  quickly  and  then  leveled  off  to  rather  constant  values  of 
Oy,  and  of  32.3  m,  38.7  m,  and  7.95  m,  respectively.  This  termination  of 
the  particle  cloud's  growth  is  due  to  the  fact  that  when  particles  follow  pure 
KS  flows  they  remain  on  particular  evolutionary  streamlines  and  are  not  able  to 
diffuse  across  to  neighboring  streamlines.  This  was  demonstrated  by  adding  a 
very  small,  high  frequency,  homogeneous  turbulence  component  having  a  velocity 
standard  deviation  of  only  o  -  4.1  cm/s  and  a  Lagrangian  time  scale  of  Tl  -  1  s. 
Such  a  level  of  turbulence  alone  would  lead  to  Oy,  o,  «  3.5  m  after  1-h,  but 
in  combination  with  the  KS  flow  fields  leads  to  standard  deviations  of  53.9  m, 
46.7  m,  and  30  m,  respectively.  Thus,  as  anticipated,  this  very  small  amount  of 
diffusion  caused  the  particle  cloud  to  grow  much  larger  than  computed  from  the 
quadrature  addition  of  uncorrelated  processes.  We  also  should  note  that  due  to 
the  fact  that  Tj,  <  At  for  the  fine  scale  turbulence  component,  the  Langevin 
equation  reduced  effectively  to  the  random  walk  problem  so  the  simpler  random 
walk  equation  was  utilized.  In  an  operational  model,  a  separate  Langevin 
equation  would  be  carried  for  each  particle  to  maintain  complete  model 
generality. 

4.3  Additional  Modules  and  Computational  Considerations 

Section  3.1  describes  the  results  of  testing  and  intercomparing  several  dry 
deposition  velocity  modules  for  particulates.  The  algorithm  extracted  from  the 
most  recent  version  of  the  Urban  Airshed  Model  (that  is,  UAM-V)  appears  to 
provide  reasonable  predictions  up  to  particle  diameters  of  about  30  /^m,  but  the 
work  of  Noll  and  Fang  (1989)  suggests  that  yet  larger  particles  (for  example,  up 
to  dp  -  100  pm)  may  deposit  with  effective  deposition  velocities  significantly 
different  from  that  estimated  by  Stokes  law. 

Section  3.2  describes  a  method  for  assigning  survival  probabilities  or, 
equivalently,  remaining  mass  fractions  to  particles  that  are  dry  deposited.  This 
approach  should  be  reasonable  for  those  particle  sizes  whose  trajectories  are  not 
significantly  altered  by  gravitational  settling.  For  larger  particles  with  more 
significant  gravitational  settling  velocities,  additional  tests  are  required  to 
ensure  that  the  survival  probabilities  computed  via  equation  (74)  do  not  cause 
any  "double  counting"  of  gravitational  settling  losses. 

A  scheme  for  computing  point  and  path-averaged  concentrations  based  on  a  modified 
kernel  formulation  is  discussed  in  section  3.3.  While  the  simple  expression 
given  by  equation  (75)  is  easily  coded,  sensitivity  studies  of  statistical 
uncertainty  versus  length  scale,  a,  and  particle  release  rate  should  be  conducted 
to  ensure  a  properly  engineered  model  with  well  understood  uncertainty 
characteristics . 

In  addition,  a  complete  multispectral  transmission/opacity  model  will  ultimately 
require  an  electromagnetic  scattering  module.  These  exist  in  a  variety  of  forms 
ranging  from  simple  optical  depth  calculations  to  full,  multiple  scattering 
computations  for  particles  having  an  arbitrary,  complex  (that  is,  real  and 
imaginary  parts)  index  of  refraction  as  a  function  of  incident  wavelength  A. 
This  topical  area  was  considered  beyond  the  intended  scope  of  the  current  Phase 
I  effort. 
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Finally,  we  note  that  the  original  proposal  suggested  that  a  RAM  disk  would  need 
to  be  utilized  to  store  the  large  number  of  particle  trajectories  anticipated  in 
an  operational  model.  The  need  arose  solely  because  of  the  640K-byte  program 
size  limit  of  the  DOS  operating  environment  on  personal  computers.  DOS  extenders 
and/or  other  operating  systems  such  as  UNIX,  are  now  more  widely  available  and 
will  enable  all  the  needed  particle  trajectory  information  to  be  retained  within 
the  program  provided  sufficient  physical  memory  {that  is,  4  to  12  mbar)  is 
available.  This  alternative  will  be  several  times  faster  than  a  RAM  disk  and 
orders  of  magnitude  faster  than  the  use  of  a  hard  disk  for  this  function. 

5.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  feasibility  of  developing  a  four - dimens ional ,  (for  example,  space  and  time 
varying)  non-Gaussian,  mesoscale  model  for  multlspectral  smoke  assessments  has 
been  evaluated.  The  basic  approach  is  not  unlike  that  suggested  by  Ohmstede  and 
Stenmark  (1980)  except  that  the  advent  of  kinematic  simulation  (KS)  now  enables 
one  to  develop  nondive r gent ,  sub -grid- scale  (SGS)  flow/ turbulence  fields  that 
evolve  in  time  so  that  true  snapshots  in  time  of  particle  distributions  may  be 
obtained  in  contrast  to  the  ensemble  average  measures  obtainable  from  use  of  the 
Langevin  equation  to  emulate  all  the  subgrid  scales  of  turbulence. 

The  basic  findings  of  this  study  include; 

(a)  Adequate  information  about  the  convective  boundary  layer  is  available 
from  analyses  of  field  data  plus  LES  data  bases  to  construct  the  appropriate  SGS 
velocity  fields  thru  kinematic  simulation. 

(b)  Available  information  on  turbulence  in  the  stable  boundary  layer  and 
that  associated  with  SGS  obstacles  is  less  complete  but  is  sufficient  for 
developing  the  SGS  velocity  field  parameters  needed  by  the  KS  approach. 

(c)  The  KS  approach  of  Fung  et  al.  (1990)  has  been  tested  for  velocity 
fields  appropriate  for  homogeneous  turbulence  in  an  unbounded  domain.  The 
methodology  achieves  the  desired  basic  objectives,  but  modified  fields  are  needed 
for  inhomogeneous  turbulence  in  an  environment  bounded  by  the  ground.  A 
plausible  set  of  KS  field  equations  have  been  developed  during  this  study  and  are 
given  by  equations  (49)  through  (51);  however,  they  have  not  been  tested. 

(d)  Several  deposition  velocity  algorithms  for  particles  have  been  tested 
and  the  resistance  based  model  contained  within  the  UAM-V  is  recommended  for 
inclusion  in  the  multispectral  smoke  model.  Additional  work  may  be  required  to 
adequately  model  deposition  velocities  for  particles  larger  than  30  /im  in 
diameter. 

(e)  A  model  for  assigning  survival  probabilities  to  individual  particles 
to  account  for  surface  removal  has  been  developed  during  this  study.  Equation 
(74)  presents  this  model  which  is  a  function  of  deposition  velocity,  time  step, 
the  harmonic  mean  height  above  terrain,  and  the  depth  of  the  lowest  model  layer 
over  which  deposition  removal  can  occur.  The  new  model  should  be  evaluated, 
particularly  for  larger  particle  sizes  where  the  particle's  trajectory  is 
significantly  influenced  by  gravitational  settling. 
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(f)  Particles  are  transported  by  a  wind  velocity  that  is  a  composite  of  the 
mesoscale  mean  flow,  the  KS  flow  associated  with  the  resolved  scales  of 
turbulence,  and  a  Langevin  equation  generated  component  associated  with  the  finer 
scales  of  turbulence  that  cannot  be  economically  accounted  for  by  the  KS 
representation . 

(g)  A  modified  "kernel"  method  is  proposed  for  ccaputing  point  and  path 
averaged  concentrations  as  a  ftinction  of  particle  size.  A  smoothing  kernel  with 
a  time- independent  length  scale  centered  on  the  sampling  point  or  line,  rather 
than  centered  on  the  particles,  is  envisioned  as  the  method  which  will  be  most 
compatible  with  the  overall  modeling  objective  of  obtaining  instantaneous 
snapshot  measures  rather  than  ensemble  or  time  average  measures.  This 
straightforward  sampling  methodology  does  not  require  significant  testing  but 
should  be  engineered  to  reflect  tolerable  levels  of  statistical  uncertainty  given 
an  operational  design  level  number  of  particle  trajectories. 

The  work  performed  during  this  Phase  I  effort  has  provided  a  number  of  discrete 
modeling  components  along  with  some  very  specific  areas,  as  mentioned  above,  for 
improving  these  components.  Nevertheless,  the  completion  of  an  operational  model 
as  depicted  in  figure  1  will  require  additional  developmental  work. 

One  major  future  work  area  involves  improvement  of  the  linkage  between  the 
current  parameterizations  of  turbulence,  as  discussed  in  sections  2.1  and  2.2., 
and  the  KS  parameters  of  section  2.4.  Fung  et  al.  (1990)  give  some  guidance  on 
the  partitioning  of  energy  into  various  spectral  components  and  consider 
plausible  spectral  shapes;  however,  the  results  of  our  section  4.2  KS  tests  show 
that  the  energy  spectrum  resulting  from  KS  modeling  is  significantly  different 
from  the  assumed  input  spectrum.  This  is  because  the  random  components  of 
amplitude  and  phase  generate  higher  and  lower  frequency  components.  In  fact,  the 
synthesized  spectra  contain  energy  at  frequencies  all  the  way  to  the  Nyquist 
cutoff.  Beyond  a  trial  and  error  approach  or  an  approach  involving  iterative 
nonlinear  optimization,  we  need  to  develop  a  clearer  understanding  of  the  linkage 
between  input  and  output  KS  spectra.  This  is  also  necessary  if  we  are  to 
understand  how  much  energy  and  what  frequency  range  should  be  accounted  for  by 
the  Langevin  particle  motion  component  of  section  2.3.  Perhaps  only  those 
frequencies  beyond  the  KS  field  advancement  frequency  need  be  represented  by  the 
Langevin  equation. 

A  Phase  II  effort  to  produce  an  operational  scheme  should  also  be  governed  by 
clear  guiding  considerations  such  as: 

•  The  completed  model  should  be  easy  to  operate  and  its  predictions 
straightforward  to  interpret.  The  operational  model  should  be  able  to 
run  on  time  scales  of  minutes  on  a  386/486  generation  personal  computer 
with  4-16  MB  of  physical  memory. 

•  The  completed  model  should  function  on  a  variety  of  computer  platforms 
with  various  operating  systems.  The  current  FORTRAN  77  programs  can 
run  under  DOS  or  UNIX;  however,  it  may  be  necessary  to  recode  the 
algorithms  into  Ada  and  perform  graphics  functions  under  XWindows  to 
achieve  maximum  compatibility  with  DoD  standards. 
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•  The  completed  model  should  produce  a  variety  of  useful  statistical  and 
graphical  outputs.  As  instantaneous  snapshot  quantities  are  computed, 
it  would  seem  reasonable  to  archive  the  entire  probability  density 
distribution  of  predictions  in  addition  to  the  time -averaged  moment 
quantities . 

Finally,  an  operational  model  should  be  designed  to  Interface  to  modules  and 
models  developed  or  used  by  the  U.S.  Army.  For  example,  the  current  prototype 
accesses  wind  and  meteorological  fields  produced  by  the  CALMET  model.  This 
interface  should  be  kept  relatively  general  so  that  meteorological  models 
possibly  better  suited  to  the  Army's  needs,  can  be  easily  accessed.  Another 
clear  example  of  desired  flexibility  involves  the  atmospheric  optics  module.  A 
multispectral  radiative  transfer  or  opacity  model  is  needed  to  convert  the 
concentration  predictions  into  atmospheric  transmissions  at  various  wavelengths. 
This  module  has  not  been  addressed  during  this  Phase  I  study  but  has  been  part 
of  other  U.S.  Army  research  programs.  Thus,  a  successful  Phase  II  model 
development/delivery  program  will  involve  a  more  detailed  evaluation  of  the 
Army's  completed  model  requirements,  an  assessment  of  the  desirability  of 
interfacing  to  available  U.S.  Army  modeling  components,  as  well  as  completion  of 
the  research  and  development  tasks  identified  above . 
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